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NOTATION 

Usage  Meaning 

•         dx 

x        — :  total  derivative  of  the  quantity  x  with  respect  to  time  t. 
dt 

x         x  is  a  column  vector. 

T 

Transpose  of  x;  defined  only  when  x  is  a  vector  or  a  matrix. 

dx 

<5— ;  partial  derivative  of  x  with  respect  to  y. 

dy 

Partial  derivative  of  the  column  vector  x  with  respect  to  the 
scalar  y.   The  result  is  a  column  vector  whose  ith  component 
is  the  partial  derivative  of  the  ith  component  of  x  with 
respect  to  y. 

dx 
x,  ,  ~=-   Partial  derivative  of  the  scalar  x  with  respect  to  the  column 

^   -   vector  y;  also  called  the  gradient  of  x  with  respect  to  y. 

It  is  a~row  vector,  whose  ith  element  is  the  partial  deriva- 
tive of  x  with  respect  to  the  ith  element  of  y. 
3x 
x,  ,  -j—   Partial  derivative  of  a  column  vector  x  by  a  column  vector  y. 
Z        ¥.        The  result  is  a  matrix  whose  (i,j)th  element  is  the  partial 
derivative  of  the  ith  element  of  x  with  respect  to  the  jtn 
element  of  y. 

T      n    2 
Norm  of  the  column  vector  x.   Defined  as  either  x  x  or   E  K.x 

~  "    i=l  1  1 

(to  be  specified  which  one)  where  x  is  a  vector  of  dimension  n, 
K.  are  given  positive  numbers,  x.  is  the  itn  element  of  x. 
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Questions  about  applicability  of  analytical  mechanics  and 
usefulness  of  optimal  control  theory  in  determining  optimal  human 
motions  arise  quite  naturally,  and  especially,  in  the  context  of  man's 
increased  activities  in  outer  space  and  under  water.   So  far  very 
little  work  has  been  done  to  answer  these  questions.   In  this  disserta- 
tion investigations  to  answer  these  questions  are  presented. 

A  particular  gymnastic  maneuver,  namely,  the  kip-up  maneuver  is 
examined  experimentally  and  theoretically.   A  mathematical  model  for 
a  human  performer  is  constructed  for  this  maneuver  from  the  best  per- 
sonalized inertia  and  joint  centers  model  of  a  human  being  available 
today.   Experiments  with  the  human  performer  and  photographic  data  col- 
lection methods  are  discussed.   Comparisons  of  the  observed  motion  with 
solutions  of  the  mathematical  model  equations  are  presented.   Discrepan- 
cies between  the  actual  motion  and  the  computed  motion  are  explained  in 
terms  of  principles  of  mechanics  and  errors  in  measurements.   Some 
changes  in  the  model  are  suggested. 


IX 


An  approximate  analytic  solution  of  the  kip-up  maneuver  performed 
in  minimum  time  is  obtained  for  the  model  via  optimal  control  theory. 
Several  numerical  methods  are  used  to  determine  the  solution,  which  is 
compared  with  the  observed  performance  of  the  human  subject.   Difficul- 
ties in  solving  human  motion  problems  by  existing  numerical  algorithms 
are  discussed  in  terms  of  fundamental  sources  of  these  difficulties. 

Finally,  recommendations  for  immediate 'future  work  have  been 
made. 


CHAPTER  1 


INTRODUCTION 


1.0.  Why  and  What 

Man' s  increased  interest  in  the  exploration  of  space  and  the 
oceans  was  an  impetus  for  a  better  understanding  of  the  mechanics  of 
large  motion  maneuvers  performed  by  human  beings.   Experience  in  space 
walking  and  certain  athletic  events  brought  out  the  fact  that  human 
intuition  does  not  always  give  correct  answers  to  questions  on  human 
motion.   For  certain  problems  the  solutions  must  be  found  by  analytical 
methods  such  as  methods  of  analytical  mechanics  and  optimal  control. 

Broadly  speaking,  this  work  deals  with  the  application  of  the 
principles  of  mechanics  and  optimal  control  theory  in  the  analytical 
determination  of  human  motion  descriptors. 

1.1.  Dynamics  and  Optimization  in  Human  Motion 

Dynamics  provides  the  basic  foundation  of  the  analytic  problem 
while  optimal  control  theory  completes  the  formulation  of  the  mathemat- 
ical problem  and  provides  means  to  solve  the  problem.   In  any  endeavor 
of  analytic  determination  of  a  human  motion  the  first  steps  are  construc- 
tion of  a  workable  model  having  the  same  dynamics  of  the  motion  as  that 
of  the  human  performer,  and  obtaining  the  equations  of  motion  for  the 
human  model.   However,  the  principles  of  mechanics  alone  do  not  give 
enough  information  for  analytical  determination  of  a  desired  maneuver 
whenever  there  are  sufficient  degrees  of  freedom  of  movement.   Without 


knowing  what  goes  on  in  the  human  motor  system,  optimal  control  is 
presently  the  only  known  analytical  method  which  can  provide  the  remain- 
ing necessary  information. 

After  a  workable  dynamic  model  of  the  human  performer  has  been 
obtained,  the  remaining  part  of  the  analytic  determination  of  a  physical 
maneuver  is  a  problem  of  control  of  a  dynamic  system  where  the  position 
vectors  and/or  orientation  of  the  various  elements  of  the  model  and 
their  rate  of  change  with  respect  to  time  (representing  the  "states"  of 
the  "dynamic  system"  being  considered)  is  to  be  determined.   The  state 
components  change  from  one  set  of  (initial)  values  to  another  set  of 
(final)  values  at  a  later  time.   The  "control  variables,"  that  is,  the 
independent  variables  whose  suitable  choice  will  bring  the  change  are 
torques  of  the  voluntary  muscle  forces  at  the  joints  of  the  various 
limbs.   However,  the  problem  formulation  is  not  mathematically  complete 
with  the  above  statements  because  there  would,  in  general,  be  more  than 
one  way  of  transferring  a  system  from  one  state  to  another  when  such 
transfers  are  possible.   Constraints  are  required  in  the  complete  for- 
mulation.  The  concept  of  optimization  of  a  certain  physically  meaning- 
ful quantity  during  the  maneuver  arises  naturally  at  this  point.   A  cost 
functional  to  be  maximized  or  minimized  gives  a  basic  and  needed  struc- 
ture to  the  scheme  for  exerting  the  control  torques.   It  may  be  expected, 
logically,  that,  unless  given  special  orders,  a  human  being  selects  its 
own  performance  criterion  for  optimization  while  doing  any  physical 
activity.   Some  of  the  physically  meaningful  quantities  that  may  be 
optimized  during  a  physical  activity  are  the  total  time  to  perform  the 
activity,  and  the  total  energy  spent  during  the  activity. 


1.2.   Previous  Work 

Very  little  work  has  been  reported  in  the  literature  so  far 
where  the  optimization  considerations  have  been  used  in  the  study  of 
human  motion.   Earlier  work  in  the  study  of  the  mechanics  of  motion  of 
living  beings  was  done  primarily  from  the  view  of  grossly  explaining 
certain  maneuvers  modeling  the  applicability  of  principles  of  rigid  body 
mechanics.   Most  of  this  work  was  done  under  either  free  fall  or  zero 
gravity  conditions.   The  models  were  made  up  of  coupled  rigid  bodies 
and  conservation  of  angular  momentum  in  the  absence  of  external  torques 
was  the  most  often  used  principle  of  mechanics. 

The  righting  maneuver  of  a  free  falling  cat  in  midair  attracted 
the  attention  of  several  authors  in  the  early  days  of  studies  of  living 
objects.   Marey's  [1]*  photographs  of  a  falling  cat  evoked  discussions 
in  1894  in  the  French  Academie  des  Sciences  on  whether  an  initial  angu- 
lar velocity  was  necessary  in  order  to  perform  the  righting  maneuver. 
Guyou  [2]  modeled  the  cat  by  two  coupled  rigid  bodies  and  explained  the 
phenomenon  with  the  aid  of  the  angular  momentum  principle  with  the  angu- 
lar momentum  of  the  entire  cat  identically  equal  to  zero.   Later,  more 
photographic  studies  were  made  by  Magnus  [3]  and  McDonald  [4,5,6]. 
McDonald  made  an  extensive  study  of  the  falling  cats  with  a  high  speed 
(1500  frames/second)  motion  picture  camera.   His  description  of  their 
motion  added  many  details  to  previous  explanations.   McDonald  found  no 


*Numbers  in  brackets  denote  reference  numbers  listed  in  the 


List  of  References. 


evidence  for  the  simple  motion  of  Magnus.   In  addition  he  studied  the 
eyes  and  the  vestibular  organ  as  motion  sensors. 

Amar  [7]  made  one  of  the  most  complete  of  the  early  studies  of 
human  motor  activities  in  1914.   This  study  of  the  relative  motion  of 
the  head,  limbs  and  major  sections  of  the  trunk  was  made  with  a  view 
to  study  the  efficiency  of  human  motion  in  connection  with  industrial 

labor. 

Fischer  [8]  considered  the  mechanics  of  a  body  made  up  of  n  links 
and  obtained  equations  of  motion  without  introducing  coordinates.   He 
made  discussions  of  applications  of  his  theory  to  models  of  the  human 
body,  but  did  not  give  applications  for  the  equations  of  motion  he  had 

obtained. 

In  recent  years  most  of  the  analytical  studies  of  human  motion 
have  been  associated  with  human  beings  in  free  fall  as  applied  to 
astronauts  maneuvering  in  space  with  and  without  external  devices. 
McDonald  [9]  made  extensive  experimental  studies  of  human  motions  such 
as  springboard  diving  and  the  "cat-drop"  maneuver.   McCrank  and  Segar 
[10]  considered  the  human  body  to  be  composed  of  nine  connected  parts. 
They  developed  a  procedure  for  the  numerical  solution  of  their  very 
complex  equation  of  motion.   Although  some  numerical  results  were  pre- 
sented, no  general  conclusions  were  drawn. 

The  most  significant  contribution  to  the  application  of  rational 
mechanics  to  problems  in  the  reorientation  of  a  human  being  without  the 
help  of  external  torques  was  made  by  Smith  and  Kane  [11].   Specifically, 
they  considered  a  man  under  free  fall.   In  this  paper  the  authors  pointed 
out  that  the  number  of  the  unknown  functions  exceeded  the  number  of  the 


equations  of  motions  that  were  obtained  for  the  system  and  recognized 
the  need  for  optimization  considerations.   In  order  to  get  the  adequate 
number  of  equations,  they  introduced  a  cost  functional  to  be  optimized 
which  consisted  of  an  integral  over  the  total  time  interval  of  some 
suitable  functional  of  the  undetermined  generalized  coordinates. 
Optimization  of  this  functional  became  a  problem  of  calculus  of  varia- 
tions, which  yielded  the  necessary  number  of  additional  equations  (the 
Euler-Lagrange  equations)  to  solve  the  original  problem  completely. 

The  approach  of  Smith  and  Kane  suffers  from  one  major  drawback- 
it  ignores  the  internal  forces  of  the  system.   The  internal  forces  due 
to  muscle  groups  at  the  various  joints  of  the  body  segments  are  mostly 
voluntary,  have  upper  bounds  in  their  magnitudes  and  are  responsible 
for  the  partially  independent  movements  of  the  various  limbs.   Because 
they  are  internal  forces,  it  is  possible  to  eliminate  them  completely 
in  any  equations  of  motion.   These  can  be  obtained,  for  example,  if  the 
entire  system  is  considered  as  a  whole.   However,  such  equations  will 
be  limited  in  number  (at  most  six,  three  from  the  consideration  of 
translation  and  three  from  the  consideration  of  rotation)  and  will,  in 
general,  be  less  than  the  number  of  the  unknown  functions.   The  intro- 
duction of  a  cost  functional  will  yield  via  the  calculus  of  variations 
the  proper  number  of  equations.   However,  the  maneuver  obtained  may  be 
beyond  the  physical  ability  of  the  individual.   It  is  therefore  essen- 
tial to  recognize  the  role  of  all  the  internal  voluntary  forces  that 
come  into  play  during  a  physical  maneuver. 

Ayoub  [12]  considered  an  optimal  performance  problem  of  the 
human  arm  transferring  a  load  from  one  point  on  a  table  to  another  point, 


The  motion  considered  was  planar.   Internal  forces  and  constraints  on 
the  stresses  were  considered.   A  two-link  model  was  considered  for  the 
arm  and  numerical  solutions  were  obtained,  using  the  methods  of  Linear 
Programming,  Geometric  Programming,  Dynamic  Programming  and  simulation. 
The  performance  criterion  was  a  mathematical  expression  for  the  total 
effort  spent  during  the  activity.   The  motion  considered  was  simple 
from  the  point  of  view  of  dynamics.   Also,  it  required  less  than  ten 
points  to  describe  the  entire  motion.   This  allowed  consideration  of 
many  physical  constraints  but  the  dynamics  of  the  human  body  motion  con- 
sidered was  quite  different  than  when  large  motion  of  the  various 
limbs  are  involved.   A  list  of  references  of  works  on  human  performance 
from  the  point  of  view  of  Industrial  Engineering  is  given  in  Ayoub's 

thesis. 

The  research  in  the  area  of  free  fall  problems  showed  that 
analytical  solutions  agreed  qualitatively  with  observations  in  those 
cases  in  which  rational  mechanics  has  been  applied  with  care.   Examples 
of  such  problems  are  the  cat-drop  problem  and  the  jack-knife-diver 
maneuver  (Smith  and  Kane).   Attempts  at  analytical  solutions  of  activ- 
ity problems  which  are  not  in  the  free  fall  category  have  not  been  com- 
pletely successful.   Such  problems  are  in  athletic  activities  (such  as 
running,  part  of  the  pole  vault  maneuver  and  part  of  the  high  jump 
maneuver)  and  in  activities  associated  with  working  on  earth. 

Simultaneously  with  the  study  of  the  dynamics  of  motion,  several 
investigations  were  made  for  the  determination  of  the  inertia  parameters 
of  human  beings  at  various  configurations.   Knowledge  of  inertia  param- 
eters is  essential  for  performing  any  dynamic  analysis.   A  list  of  the 


research  activities  in  this  area  is  given  in  References  13,  14  and  15. 
However,  only  Hanavan  [14,15]  has  proposed  a  personalized  model  of  a 
human  being.   This  inertia  model  has  been  used  in  the  present  investi- 
gation. 

1.3.   The  Problem  Statement 

The  present  work  belongs  to  a  broader  program  whose  objective 
is  to  investigate  the  basic  aspects  of  the  applicability  of  rational 
mechanics  to  the  solutions  of  any  human  activity  problem.   The  aspects 
presently  being  considered  are: 

1.  Construction  of  an  appropriate  presonalized  model  for  the 
individual's  maneuver  under  consideration 

2.  Formulation  of  a  well-posed  mathematical  problem  for  the 
analytic  description  of  the  maneuver 

3.  Solution  of  the  mathematical  problem  .by  a  suitable  analytical 
method 

4.  Comparison  of  the  analytical  solution  with  an  actual  motion 
conducted  in  an  experiment  with  a  human  subject 

5.  Determination  of  muscle  activity  and  comparison  with  computed 
muscle  torque  histories  for  the  maneuver. 

If  a  complete  analysis  based  on  the  above  steps  results  in  the 
correct  motion  as  compared  with  the  experiment,  then  the  results  can  be 
used  for  training  purposes  and  design  of  man-machine  systems,  but,  more 
importantly,  the  results  will  establish  the  applicability  of  rational 
mechanics  to  the  solution  of  problems  of  human  activity. 

In  the  present  investigation,  a  particular  gymnastic  maneuver, 
the  kip-up,  has  been  selected  for  analysis  as  outlined  above.   The  methods 
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developed  will,  of  course,  be  valid  for  other  maneuvers  but  the  ana- 
lytical model  will  be  personalized  to  the  maneuver  and  individual. 
Among  all  the  common  physical  and  athletic  activities,  the  kip-up 
maneuver  was  found  to  be  particularly  well  suited  for  analysis.   The 
motion  involves  large  motions,  and  is  continuous  and  smooth.   Also,  it 
is  planar  and  needs  relatively  fewer  generalized  co-ordinates  for  its 
complete  description,  since  a  correctly  executed  kip-up  exhibits  three 
"rigid"  links.   At  the  same  time  it  is  not  a  trivial  problem  from  the 
point  of  view  of  our  basic  objective. 

The  physical  quantity  (the  performance  criterion)  chosen  for 
optimization  (minimization,  in  this  case)  is  the  total  time  to  do  the 
maneuver. 

1.4.   The  Kip-Up  Maneuver 

The  kip-up  maneuver  is  an  exercise  that  a  gymnast  performs  on 
a  horizontal  bar.   The  gymnast  starts  from  a  position  hanging  vertically 
down  from  the  horizontal  bar  and  rises  to  the  top  of  the  bar  by  swinging 
his  arms  and  legs  in  a  proper  sequence.   During  the  maneuver  the  moti< 
is  symmetric  and  does  not  involve  bending  of  the  elbows  and  the  knees. 
Normally,  the  grip  on  the  bar  is  loose  most  of  the  time.   For  an  inex- 
perienced person,  the  maneuver  is  not  easy  to  perform. 


Lon 


1.5.   Present  Work 

In  Chapter  2,  a  mathematical  model  for  the  kip-up  motion  and 
results  of  laboratory  experiments  to  test  the  accuracy  of  the  model 
are  presented.   First,  a  mathematical  model  of  a  professional  gymnast 
for  the  kip-up  motion  has  been  constructed.   The  dynamic  equations  of 


motion  for  the  mathematical  model  were  then  obtained.   Two  sets  of 
equations  of  motion  were  derived:   one  for  the  purpose  of  verifying  the 
accuracy  of  the  mathematical  model  from  experiments  and  one  for  the 
purpose  of  optimization  of  the  kip-up  maneuver.   Next,  the  results  of 
laboratory  experiments  with  the  gymnast  are  presented.   The  gymnast  was 
told  to  perform  symmetric  maneuvers  on  the  horizontal  bar,  including  the 
kip-up.   The  maneuvers  were  photographically  recorded.   Two  of  the  many 
records  were  selected,  one  of  a  simple  swinging  motion  with  relatively 
small  oscillations,  and  another  of  his  quickest  kip-up  maneuver. 
The  angle  measurements  from  these  film  records  were  then  used  in  the 
equations  of  motion  to  check  the  accuracy  of  the  mathematical  model. 
An  error  analysis  was  then  performed  to  explain  the  disagreement  between 
the  experimental  and  the  computed  results. 

In  Chapter  3,  an  analytic  solution  of  the  minimum  time. kip-up 
for  the  mathematical  model  was  obtained  by  numerical  computations. 
First,  the  analytical  problem  of  the  determination  of  the  kip-up  in  mini- 
mum time  wa  s  stated  in  precise  mathematical  terms.   This  involved  repre- 
senting the  equations  of  motion  in  the  state  variable  form,  specifying 
the  boundary  conditions  on  the  state  variables,  establishing  the  bounds 
on  the  control  variables,  and  modeling  the  stiffness  of  the  shoulder 
and  the  hip  joints  at  extreme  arm  and  leg  movements.  A  survey  of  the 
necessary  conditions  for  time  optimality  given  by  the  optimal  control 
theory  has  been  presented.   Finally,  several  numerical  schemes  used  for 
solving  the  kip-up  maneuver  problem  are  presented  and  the  results  of  the 
numerical  computations  are  discussed. 

In  Chapter  4,  final  conclusions  of  the  present  investigation 
and  recommendations  for  future  work  have  been  presented. 


CHAPTER  2 


EXPERIMENTATION  AND  CONSTRUCTION 
OF  THE  MATHEMATICAL  MODEL 


2.0.  Introduction 

In  this  chapter,  modeling  of  a  human  being  for  the  kip-up 
maneuver  is  considered.   In  Section  2.1,  a  mathematical  model  for  a 
professional  gymnast  is  constructed,  using  the  personalized  model  of 
Hanavan  [14].   In  Section  2.2,  equations  of  motion  for  the  mathematical 
model  are  derived  for  the  purpose  of  using  them  in  the  analytic  deter- 
mination of  the  subject's  optimal  kip-up  motion.   An  equivalent  system 
of  first-order  differential  equations  are  derived  in  Section  2.3  for 
testing  the  inertia  properties  and  structure  of  the  model.   In  Sec- 
tion 2.4,  the  laboratory  experiments  performed  are  described.   The 
results  of  the  experiments  are  discussed  in  Section  2.5.   An  analysis 
of  the  comparisons  between  the  observed  and  computed  results  are  also 
presented. 

2.1.  Mathematical  Model  of  the  Kip-Up 

The  equations  of  motion  of  a  deformable  body  such  as  the  human 
body  are  usually  partial  differential  equations.   Presently,  not  enough 
is  known  or  measurable  about  the  deformation  of  the  human  body  under 
voluntary  motion  to  determine  a  partial  differential  equation  model. 
Also,  such  models  are  quite  difficult  to  handle.   However,  for  situations 
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where  the  deformation  is  small  compared  to  the  displacement  of  a  body, 
the  deformable  body  may  be  considered  as  a  rigid  body  in  writing  the 
equations  of  motion.   The  rigid  link  assumption  has  been  used  widely 
in  modeling  human  beings.   The  personalized  model  of  Hanavan  [14]  is 
based  on  this  assumption.   The  elements  of  the  Hanavan  model  are  shown 
in  Figure  1  and  consist  of  fifteen  simple  homogeneous  geometric  solids. 
This  construction  allows  a  large  number  of  degrees  of  freedom  for  the 
model  and  minimizes  the  deformation  of  the  elements  without  undue  com- 
plexity.  As  the  number  of  degrees  of  freedom  increases  to  the  maximum, 
it  is  likely  that  a  mathematical  model  based  on  the  Hanavan  model 
becomes  more  accurate.   However,  increases  in  the  degrees  of  freedom 
increase  the  model's  complexity  and  make  it  difficult  to  analyze  mathe- 
matically. 

Observation  of  a  correctly  performed  kip-up  indicates  that  the 
human  performer  might  be  modeled  quite  accurately  as  a  system  of  only 
three  rigid  links.   The  two  arms  form  one  link,  the  head-neck-torso 
system  forms  another,  and  the  two  legs  form  the  third  link.   The 
shoulder  and  hip  joints  may  be  approximated  as  smooth  hinges  where 
these  links  are  joined  together.   Deformation  of  the  link  consisting 
of  the  head,  neck,  and  torso  during  certain  periods  of  the  maneuver 
is  detectable  when  observed  via  high  speed  filming.   It  was  felt  that 
the  effect  of  the  deformation  of  the  torso  would  be  no  more  significant 
than  the  standard  deviation  in  the  Hanavan  inertia  parameters  because 
this  part  of  the  orientation  of  the  torso  is  determined  by  the  line 
joining  the  shoulder  and  the  hip  joint  centers. 
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Figure  1.   Hanavan' s  Mathematical  Model 
of  a  Human  Being. 
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The  muscle  forces  acting  at  the  hip  and  shoulder  to  cause  link 
motion  have  been  replaced  with  their  rigid  body  equivalent  resultant 
forces  and  couples.   If  the  masses  of  the  muscles  causing  the  forces 
are  small  compared  to  the  other  portions  of  the  links,  then  the  net 
effects  of  the  forces  are  the  torques  of  the  couples.   The  resultants 
do  not  appear  in  the  equations  of  motion. 

The  three-link  kip-up  model  is  shown  in  Figure  2.   It  is  con- 
structed with  elements  of  the  Hanavan  model.   Twenty-five  anthropometric 
dimensions  of  the  gymnast  were  taken  and  used  in  a  computer  program  for 
calculating  first  the  inertia  properties  of  the  Hanavan  model  and  then 
those  of  the  kip-up  model.   The  determination  of  the  inertia  properties 
of  the  kip-up  model  from  those  of  the  Hanavan  model  is  presented  in 
Appendix  A. 

2.2.   The  Equations  of  Motion 

The  mathematical  model  for  the  kip-up  motion  is  a  three-link 
system  executing  plane  motion  under  gravity.   The  active  forces  in  the 
system  are  the  pull  of  gravity   acting  on  each  of  the  elements  of 
the  system  and  the  two  muscle  torques.   The  muscle  forces  at  the  shoulder 
acting  at  the  joint  between  links  1  and  2  are  replaced  by  the  torque  u  . 
Likewise,  u   is  the  torque  for  the  hip  joint  between  links  2  and  3. 
The  system  is  suspended  from  a  hinge  at  the  upper  end  of  link  1,  repre- 
senting the  fists  gripping  the  horizontal  bar  which  is  free  to  rotate 
on  its  spherical  bearings.   All  joint  hinges  are  assumed  to  be  friction- 
less.   The  general  three-link  system  for  which  the  equations  of  motion 
will  be  derived  is  shown  with  nomenclature  in  Figure  3.   Equations  of 
motion  of  the  system  were  obtained  via  Lagrange's  equations  as  follows: 
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CGI 
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element  3 


=  length  of  element  1  =  distance  between  the  hinges  0  and  A 

=  length  of  element  2  =  distance  between  the  hinges  A  and  B 

=  distance  between  the  center  of  gravity  of  element  1  and  the  hinge  0 

=  distance  between  the  center  of  gravity  of  element  2  and  the  hinge  A 

=  distance  between  the  center  of  gravity  of  element  3  and  the  hinge  B 

,  CG2,  CG3  =  locations  of  centers  of  gravity  of  elements  1,  2,  and  3, 
respectively 

=  angle  between  O-CGl  and  OA 

=  angle  between  A-CG2  and  AB 

,  m  ,  m  =  mass  of  elements  1,  2,  and  3,  respectively 

1,1=  moments  of  inertia  of  elements  1,  2,  and  3,  respectively, 
abou?  axes  perpendicular  to  the  xz  plane  through  their  respective 
centers  of  gravity 

=  moment  of  inertia  of  the  horizontal  bar  at  the  hinge  0  about  its 
longitudinal  axis 


g  s=  acceleration  due  to  gravity 


Figure  3.   The  Three-Link  System 
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If  we  define 


12  2  ,2  2  .2  ,2 

A  =  j(  Vl  +  Ix  +  m2r2  +  W!  +  I3  +  Vs  +  Vl  +  ^2  +  V 


m2Vl 


C  =   m3l±l2 


D  =  m3r3X1 


E  =  m3r3^ 


F  =  |(V2+I2+Vl+m3r3+I3) 

J  =  m3r3+I3 


(2.2.1) 


N  =  (m2  +  m3)  i^g 


V  =  m3l2g 


W  =  m2r2g 


R  =  m3r3g 


with  Equations  (2.2.1),  we  can  express  the  Lagrangian  of  the  system  as: 
L  =  <p2[A+  B  cos  (9+3)  +  C  cos  9  +  D  cos  (9+^)  +  E  cos  f]  +  6  [F+E  cos  i|t] 

1    *  2 

+    9  cp  [2F+  B   cos    (6+P)   +    C  cos    9+  D  cos    (9+\|f)   +    2E   cos   i|r]   +    2   i|r  J 

•   •  •  * 

+    9  i|f  [ J  +  E  cos   ijr]   +  cp  i|r  [J  +  D  cos    ( 9+  \|r)   +   E  cos   \|r]   +   M  cos    (cp+aO 

+  N  cos  cp  +  V  cos  (cp+9)  +  W  cos  (cp+9+g)  +  R  cos  (cp+9+i|0-      (2.2.2) 

For  the  Hanavan  model  of  *ur  system,  the  angles  a  and  0  are  zero. 

If  we  define  un  and  u   as  positive  when  they  tend  to  increase 
the  angles  9  and  \|r  respectively,  we  can  write  the  equations  of  motion 


as: 
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_!  *L   -    *L  =    0  (2.2.3) 

dt  aj    ^~ 

d     9L   _    BL  (2.2.4) 


(2.2.5) 


dt  Be 

"  "56  = 

ui 

d     OL 

dt      • 

5l 
"  3ijF  = 

U2 

Let  us   define 

a     =  2 (A  +    (BfC)    cos    9  +   D  cos    (9+\|0    +   E   cos   l|0 

a     =  2F  +    (BfC)    cos    9  +   D   cos    (9+i|0   +    2E   cos    l|l 

a     =  J  +   D  cos    (9+t)    +   E   cos    i|i 

bl   =  a2 

b     =  2(F  +   E   cos   'jr) 

0 


b3 

S= 

J     + 

E 

cos 

Cl 

= 

a3 

C2 

= 

b3 

C3 

= 

J 

d     =  cp9[2(B+C)    sin   9  +    2D  sin    (9+t)  ]  +  9tC  2D  sin   (6+l|0   +    2E   sin  f] 
+    9t[2D  sin    (9+|)  +  2E   sin  i|t]       +       6    [(BfC)    sin   9  +    D   sin    (9+t)] 
+    ^f2[D  sin   (9+f )  +  E   sin  \|r]  -  (MfN)    sin  9  -  (V+W)    sin   (cp+9) 

-  R  sin   (<p+6+i|0 
d     =u+i[2Cp+29+i]E   sin   i{f  -cp    [(BfC)    sin   9  +   D  sin   (9+\|f)] 

-  (v+W)    sin   (9+ 9)    -  R  sin   (cp+8+\|0 

•2  *2 

d     =   u     -cp    [D  sin    (8+ijr)   +   E   sin  i|r]    -    9     E   sin  i|r  -29?      E   sin  i|r 

-  R  sin    (cp+e+ilO  . 

(2.2.6) 
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With  Equations  (2.2.6)  we  obtain  from  Equations  (2.2.2)  -  (2.2.5)  the 
equations  of  motion 


a^>  +   b16  +  c^  =  d1 


*2*  +   b26  +  C2*  =  d2 


a3cp  +  b36  +  c3t  =  d3 


(2.2.7) 
(2.2.8) 
(2.2.9) 


It  will  be  helpful  to  express  the  equations  of  motion  in 
normal  form  in  formulating  optimal  control  problems  for  the  system. 
For  that  purpose  we  define  the  state  variables  as 


X,  =  cp   X.  =  9   X, 


X4=9    X5=*    V* 


(2.2.10) 


31   bl   Cl 


a2   b2   C2 


a3   b3   C3 


(2.2.11) 


Al  = 


dl  bl  Cl 
d2  b2  C2 
d3   b3   °3 


(2.2.12) 


al   dl   Cl 


32   d2   C2 


33   d3   C3 


(2.2.13) 


31   bl   dl 


32   b2   d2 


33   b3   d3 


(2.2.14) 
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Using  these  definitions,  the  equations  of  motion  (2.2.7) 
then  be  written  in  the  normal  form  as 

*1  *  X2  -    *■ 

A, 


(2.2.9)  can 


X4  =  X 
4    A 


W 


X5=X6 


X3=X4 


6  "  A 


To  write  down  the  equations  of  motion  in  more  convenient  forms,  we  shall 
further  define  the  following  quantities: 


Al      A2      L3   T 
f(X,u)  =  [X2,  T,  X4,  — ,  X6,  T] 


ei  =  dl 


e2-d2 


Ul 


63  "  d3  "  U2 


J 


Gl 

bl 

Cl 

*1  = 

62 

b2 

C2 

63 

b3 

C3 

'*! 

Gl 

Cl 

*2  = 

a2 

G2 

C2 

31 

63 

C3 

\ 

bl 

el 

*V 

32 

b2 

C2 

a3 

b3 

°3 

(2.2.15) 


(2.2.16) 


(2.2.17) 


(2.2.18) 


(2.2.19) 
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A0   lT 


r         Ll  L2  a3  1 

a«  =  |_V  T'    x4'  T»    V  Tj 


B(X)    =   v 


'-Vs+W 


("alC3-a3Cl) 


<-alb3+a3bl) 


'Vq-W 


("alC2+a2Cl) 


'■iVaV 


>  (2.2.20) 


-/ 


V*>  =  A 


("blC3+b3Cl) 


(alC3-a3Cl} 


("alb3+a3bl) 


(2.2.21) 


V*>=a 


(blC2-b2Cl) 


(-aiC2+Vl) 


(aib2  "  a2bl} 


(2.22) 


Using   the  def initions(2. 2. 16)    -    (2.2.22),    the   equations   of  motion 
(2.2.15)    may  now  be  expressed  by   any  one  of   the  following  equations: 
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55-  "  'S^iV3  .♦s,?)-n 


5,  V,     +  S-,-?   <\-\    i'/,dA 

3.1    2,;       Jo         J.J  J 


S' t! 


J':i 


i'       t:       : 


-S;3r3rY.di  ♦  r{S3*-:s,>5V3t?+5.t3J,n,  -   ^aVj^Sj^ngjrtyA 


V  I    L  . 


J  J 


3  c   3  J      3 


*     jiS-  V-       V-,:        dV 


and 


(CI  -    f:v_rv;  «., 


V  I  t  I 


dA  +■   f{v,-V-]  ■-•$- 


v(t) 


(9) 


L'p  to  this  point.  this  analysis  has  beer.  »erv  gen- 
era" and  tr  >-es-d  iien  s  iona",  •witn  na  reference  to  a  ?~ate- 
like  two-dimensional  problem.   Equation  i"8)  «ntn  the  aid 
of  (?;  car.  be  applied  to  any  ihree-fiisensional  elasticity 
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with  cp  expressed  in  terms  of  p  using  Equation  (2.3.1).   Hamilton's 
canonic  equations  are  then  given  by 

<p  =  ||  (P,«P,t)  (2-3-3> 

and 

9H 
p  =  "  Sp 

From  Equations  (2.3.2)  -  (2.3.4)  we  obtain 


(P.tp.t).  (2.3.4) 


p-9[2F+(B+C)cos    6+D  cos  (9+'>)+2E  cos  \|f]  -  jt[J+D  cos  (9+jf)  +  E  cos  i|f] 
cp   = — — —  _— 

2[A  +    (B+C)    cos    9  +   D  cos    (8+f)    +  E   cos   i|f] 

(2.3.5) 

p  _  _[(M+N)  sin  cp  +  (V+W)  sin  (9+9)  +R  sin  (9+6+i|()].  (2.3.6) 

If  we  wish  to  introduce  the  effects  of  the  friction  at  the  hinge  0, 
the  equation  for  p  becomes 

*     *  -  f  (2.3.7) 

H     09    r 

where  F   is  the  generalized  friction  torque  at  the  hinge  0. 
r 

The  Integration  Scheme 

The  integration  scheme  integrates  Equations  (2.3.5)  and  (2.3.7) 
The  initial  condition  of  9  is  obtained  from  the  measured  values  of  9. 
The  initial  value  of  p  is  obtained  from  the  definition  of  p  given  in 
Equation  (2.3.1).  To  compute  the  initial  value  of  p ,  9  has  been  com- 
puted for  this  starting  point  only.  The  9  and  1/  data  are  differenti- 
al! .  OH 
ated  numerically  to  generate  9  and  i|i  values  at  every  step.   -^  and  -^ 

are  computed  at  every  step  by  using  these  values.   The  value  of  Fr  is 
not  known  a  priori  and  it  requires  a  separate  experiment  for  its 
determination.   Integrations  were  done  with  Fr  =  0  and  F^ =  C  sgn  9 
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with  C  determined  experimentally.   The  difference  between  the  two 
cases  was  found  to  be  insignificant  even  with  C=l  ft-lb  which  was 
well  above  the  possible  friction  torque  at  the  bearings. 

2.4.   Experimental  Procedure 

The  test  rig  consisted  of  a  horizontal  bar  and  two  motion 
picture  cameras.   The  horizontal  bar  was  made  of  a  short  solid  round 
steel  bar  1-3/16  inches  in  diameter  and  58  inches  long  supported  by  two 
very  rigid  vertical  columns  through  a  pair  of  self-aligning  spherical 
bearings.   The  bearings  allowed  free  rotation  of  the  bar  with  the  arm 
of  the  subject. 

One  movie  camera  was  placed  with  its  line  of  sight  aligning 
with  the  horizontal  bar  and  about  30  feet  away  from  the  bar  as  shown 
in  Figure  4.   Alignment  of  the  camera's  line  of  sight  with  the  horizontal 
bar  would  give  the  correct  vertical  projection  of  the  two  arms  on  the 
film  for  determining  the  angle  cp  directly.   The  second  camera  was  placed 
in  front  of  the  horizontal  bar  at  the  same  elevation  as  the  bar.   The 
film  taken  in  this  camera  showed  whether  a  particular  motion  was  sym- 
metric or  not  about  the  vertical  plane  of  motion  and  also,  the  angle 
between  the  two  arms,  which  is  required  for  computing  the  moment  of 
inertia  of  the  arms.   The  film  speed  was  determined  from  the  flashes  of 
a  strobe  light  regulated  by  a  square  wave  generator. 

Experiments 

The  experiments  were  done  during  two  separate  periods  with  the 
same  subject,  a  professional  gymnast.   An  average  of  15  experiments 
of  the  subject's  performance  on  the  bar  were  recorded  on  film  each  day. 
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The  subject  was  told  to  avoid  bending  his  arms  and  legs  and  to  maintain 
symmetric  motion.   In  the  early  days  of  experimentation  he  was  told  to 
just  swing  on  the  bar  by  moving  his  stiff  arms  and  legs  relative  to  the 
torso.   These  experiments  were  done  with  the  idea  of  obtaining  small 
angle  data  for  verifying  the  inertia  properties  of  the  Hanavan  model. 
The  filming  of  the  camera  was  done  at  speeds  of  32  frames  per  second  and 
62.1   frames  per  second  (as  determined  from  stroboscopic  measurements). 
In  later  experiments  the  subject  was  told  to  perform  the  maneuver  with 
some  specific  objectives.   He  was  told  to  perform  what  he  thought  would 
be  the  kip-up  (1)  in  minimum  time,  (2)  with  minimum  expenditure  of 
energy,  and  (3)  putting  "least  effort."  Each  of  these  maneuvers  was 
repeated  several  times.   Between  any  two  subsequent  maneuvers ,  the 
subject  was  given  adequate  rest  periods  to  avoid  fatigue.   This  exper- 
iment was  conducted  with  the  idea  of  making  optimization  studies  as  well 
as  testing  the  model. 

Four  white  tapes  were  stuck  to  the  subject,  on  the  sides  of 
his  upper  and  lower  arms,  sides  of  his  torso,  and  on  the  sides  of  his 
legs.   These  tapes  were  aligned  between  joint  centers  as  suggested  by 
the  Hanavan  inertia  model. 

Processing  the  Data 

The  film  speed  was  measured  with  the  aid  of  a  stroboscope  by 
running  the  camera  with  a  developed  film  with  the  same  speed  setting. 
After  removing  the  lens  the  shutter  was  exposed  to  the  stroboscope 
flash.   By  arresting  the  shutter  in  the  stroboscope  light,  the  shutter 
speed  was  obtained. 


26 

The  films  were  run  on  an  L-W  Photo  Optical  Data  Analyzer. 
For  the  purpose  of  testing  the  inertia  properties,  two  maneuvers  were 
selected,  one  from  each  day's  filmings.   The  films  were  projected 
plumb  line  perpendicularly  on  paper  fixed  to  a  vertical  wall.   As  each 
frame  was  projected,  the  white  tapes  fixed  on  the  subject,  now  clearly 
visible  in  the  image,  were  marked  out  on  the  paper  of  the  pad  by  means 
of  a  pencil  and  a  straight  edge.   In  this  way  each  frame  was  "trans- 
ferred" on  separate  sheets  of  paper.   The  angles  9  and  lji  were  then 
measured  from  these  traces.   The  angle  cp  was  measured  with  respect  to 
a  vertical  reference.   The  vertical  reference  was  obtained  from  a  sharp 
window  wall  in  the  background. 

Of  the  two  sets  of  data  processed,  one  was  smoothed  before 
using  it  in  the  integration  scheme.   This  was  the  one  filmed  at  the 
higher  speed  for  the  fast  kip-up  motion.   A  plot  of  the  raw  data  and 
those  after  preliminary  smoothing  for  this  set  are  shown  in  Figure  5. 

2. 5  Results  and  Discussion 

The  results  of  the  integration  of  the  equations  of  motion  for 
two  of  the  data  sets  analyzed  are  shown  in  Figures  6  and  7.   Figure  6 
shows  data  for  swinging  motion,  while  Figure  7  is  for  the  kip-up. 
Also,  in  Figure  7  is  given  the  computed  results  corresponding  to 
reinitiating  the  integration  program  with  the  measured  data.   Differ- 
ent starting  points  of  integration  were  selected  to  eliminate  the 
errors  that  were  generated  before  these  points. 

In  Figure  6,  curve  2  of  the  computed  values  for  the  unsmoothed 
data  agrees  well  with  curve  1  of  the  measured  values  for  a  little  more 
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The  results  obtained  so  far  art  applicable  to 
plates  cons;stir.c  of  any  nuiitt-e*-  or'  layers  having  differ- 
ent properties  an;  geometry.   F!jt  "-on.  '.ere  en,  the  analy- 
sis would  depend  -n  the  type  of  plate  considered  vi  z .    sin- 
gle layer  or  mu'ti player  elates,  thin  or  moderately  thick 
plates,  etc.,  and  also  on  the  type  of  displacement  func- 
tions V  and  V,   Habip  [41]  has  demonstrated  the  applica- 
x  3  " 

oility  of  these  results  to  a  single  layer  plate, 

in  this  dissertation,  a  three-layered  plate,  popu- 
larly known  as  a  sandwich  plate  will  be  considered. 

Sandwi  ch  ^1  ate 


A  sandwich  plats  consisting  of  three  layers  is 
shown  in  Figure  (2).   T"e  face  layers  are  much  thinner 
th  =  r.  the  core.   All  layers  are  uniformly  thick  throughout- 
The  two  faces  are  of  the  sacie  thickness  t<,  the  core  thick- 
ness being  t  . 

Tie  present  analysis  is  capable  of  creatine  nixed 
boundary  value  problems.   On  the  upper  surface,  the  line 
AB  separates  the  two  regions  -Avf  -  and  ,A  ,y,  over  which 
displacements  and  stresses  respectively  are  prescribed. 
On  toe  lower  surface,  the  line  CD  separates  the  correspond- 
ing regions.   'he  line  CD  is  located  exactly  below  the  line 
A.B.   Also,  at  any  two  points  located  or  the  up  sen  and  lower 
surfaces  of  the  plate  and  having  t n e  s a m e  X ,  and  x  - 
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2.6.   Sources  of  Errors 

The  following  are  considered  to  be  the  sources  of  errors 
responsible  for  the  disagreement  between  the  computed  and  measured 
values. 

2.6.1.  Imperfections  in  the  Model 

The  human  being  for  the  motion  studies  was  modeled  as  a  system 
of  rigid  bodies.   The  response  of  the  system  to  be  compared  to  that  of 
the  rigid  body  model  was  a  single  generalized  coordinate  of  the  system. 
The  errors  in  modeling  can  be  lumped  into  the  overlapping  categories  of 
(1)  definition  of  the  generalized  coordinates  of  the  rigid  elements  of 
the  system,  (2)  deformations  of  link  lines  during  motions  from  prior 
joint  center  measurements,  and  (3)  significant  variations  during  motion 
in  the  inertia  properties  of  the  torso  with  respect  to  the  fixed  coor- 
dinate system.   In  several  experiments  the  torso  deformed  with  signif- 
icant movement  of  the  shoulder  joint  centers.   In  these  cases  the  con- 
stant inertia  properties  model  is  obviously  incorrect.   These  variations 
not  only  cause  inaccuracies  in  the  inertia  parameters  but  also  result  in 
errors  in  the  link  lines  or  additional  errors  in  the  mass  center  of 
element  2.   These  errors  are  reflected  in  errors  in  the  angles  8  and  i|r, 
which  in  turn  cause  a  time  varying  "phase  shift"  in  the  computed  angle  cp. 

2.6.2.  Errors  in  Filming  and  Processing  the  Data 

The  primary  source  of  errors  in  filming  was  considered  to  be 
caused  by  inaccuracies  in  knowing  precisely  where  the  link  lines  were 
in  relation  to  the  film  plane.   As  mentioned  in  Section  2.4,  care  was 
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taken  to  minimize  this  error.   That  is,  the  cameras  were  aligned  with 
the  horizontal  bar  so  that  the  link  line  of  the  arms  projected  in  the 
plumb  line  plane  of  the  film  was  very  nearly  the  correct  model  refer- 
ence line.   Since  the  link  lines  of  the  torso  and  legs  were  nearly 
plumb  line  vertical,  no  corrections  of  the  image  data  were  required 
for  these  links.   Filming  of  static  thin  rods  which  were  connected  to 
the  bar  at  known  angles  to  each  other  and  the  vertical  plane  produces 
overall  measurement  errors  of  about  1  degree  standard  deviation. 
Errors  up  to  3  degrees  can  be  expected  in  data  from  films  of  the  motion 
experiments.   These  errors  can  cause  the  rates  of  the  angles  to  be  in 
error  by  more  than  30  percent.   This  is  the  primary  cause  of  the  error 
in  the  amplitude  of  the  angle  9  computed  from  the  dynamical  equations. 
(Second  derivatives  of  the  data  can  be  in  error  by  more  than  100  percent. 
This  ruled  out  the  use  of  equations  such  as  Euler's  or  Lagrange's.) 

2.6.3.   The  Integration  Scheme 

The  integration  scheme  uses  at  any  step  the  measured  values  of 
0  and  i|;  and  the  computed  values  of  9  and  i|r  stored  previously.   Once  a 
difference  between  the  measured  (actual)  and  the  calculated  values  of  cp 
has  developed  at  a  time  t  due  to  any  of  the  sources  of  errors  discussed 
above,  the  system  configuration  determined  by  the  calculated  value  of  cp 
and  the  measured  (actual)  values  of  9  and  t|r  at  the  time  t  will  be  dif- 
ferent from  that  of  the  actual  system  at  that  time.   This  will  cause 
the  model  to  have  a  different  response  after  time  t  than  that  of  the 
actual  system  which  has  a  different  relative  configuration.   This 
in  turn  will  cause  a  further  deviation  between  the  actual  and  the 
calculated  motion  that  follows  this  instant  of  time.   To  reduce  this 
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effect  of  propagation  of  error  via  the  system  equations,  several 
reinitializations  of  cp  and  p  were  done  by  restarting  the  integrations 
at  different  points. 

Various  order  differentiation  and  integraion  schemes  were 
tested  with  insignificant  differences  in  the  results  for  smoothed 
data.   These  were  done  with  data  from  filming  at  F=62.1  frames  per 
second  with  a  maximum  integration  step  size  of  2/F  second.   As  men- 
tioned previously,  the  main  cause  of  amplitude  error  was  the  errors  in 
the  derivatives  of  the  raw,  unsmoothed  data. 

The  results  obtained  from  the  experiments  show  that  the  model 
for  the  kip-up  motion  constructed  from  the  Hanavan  model  was  reason- 
ably good  considering  its  kinematical  simplicity.  In  spite  of  imper- 
fections in  the  model,  its  dynamic  behavior  was  quite  similar  to  the 
actual  motion,  so  that  this  model  could  provide  reasonable  estimates 
of  optimal  human  performance  via  the  theory  of  optimal  processes  and 
numerical  solution  methods. 

However,  the  results  obtained  from  the  experiments  indicate 
that  the  application  of  rational  mechanics  to  the  analysis  and  design 
of  man-machine  systems  could  prove  inadequate  unless  the  model  and  the 
data  gathering  techniques  can  be  improved.   This  is  especially  true  in 
the  design  of  high  accuracy  or  low  tolerance  systems. 


CHAPTER  3 


ANALYTIC  DETERMINATION  OF 
THE  MINIMUM -TIME  KIP-UP  STRATEGY 


3.0.   Introduction 

In  this  chapter  the  determination  of  an  analytic  solution  of  the 
kip-up  maneuver  is  presented.   The  problem  of  analytical  determination 
of  the  kip-up  strategy  in  minimum  time  has  been  cast  as  a  problem  of 
optimal  control  of  dynamical  systems.   Before  the  techniques  of  the 
optimal  control  theory  may  be  applied  to  the  problem,  it  is  necessary 
to  state  the  physical  problem  in  the  language  of  mathematics  and  to 
introduce  the  physical  constraints  that  must  also  be  considered  for  the 
solution.   Thus,  the  first  four  sections  of  this  chapter  have  been 
devoted  to  the  formulation  of  the  mathematical  problem.   In  Section  3.5 
a  survey  of  the  necessary  conditions  for  optimal ity  obtained  from  the 
optimal  control  theory  is  presented.   Since  the  problem  under  consider- 
ation cannot  be  solved  in  closed  form,  numerical  methods  were  used  to 
obtain  the  solution.   In  Sections  3.6  -  3.9,  the  choice  of  the  numerical 
methods,  their  derivations  and  the  results  of  the  numerical  computations 
are  discussed.   In  Section  3.10,  results  of  the  numerical  computations 
are  compared  with  the  actual  motion. 
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3.1.   Mathematical  Formulation  of  the  Kip-Up  Problem 

The  problem  is  to  determine  the  minimum  time  strategy  for  the 
man  model  to  kip-up  without  violating  control  constraints.   These  con- 
straints represent  the  maximum  torques  the  man's  muscles  can  exert  for 
any  given  configuration.   Formulated  mathematically,  we  have  the  follow- 
ing: 

For  the  system  equations 

£  =  t(h^)    =  £(P  +  B<29H  (3.1.1) 

and  the  boundary  conditions 

X(0)  =  X       (given)  (3.1.2) 

—      —  o 

§(X(tf))^  X(tf)  -  Xf  =  0  ,  xf  =  given  ■  (3>1>3) 


where 


and 


given  by 


t  =  final  time,  to  be  determined  (3.1.4) 


X(t)  is  the  time-dependent  state  vector, 


xx(t)  =  cp(t),   x2(t)  =  cp(t),   x3(t)  =  9(t) 


X„(t)  =   9(t),   X_(t)  =   ^(t),   X  (t)  =  i|r(t) 
4  5  t> 


(3.1.5) 


find  a  control 

u(t)  =  [u1(t),u2(t)]T  .    (3.1.6) 

such  that  simultaneously  Equations  (3.1.1)  -  (3.1.3)  are  satisfied, 

t   is  minimized,  and  for  all  values  of  t,  0  S  t  S  tf)  the  inequalities 


S*(X)  £  ux(t)  *  S^(X) 


S,(X)  £  u  (t)  £  S^(X) 


(3.1.7) 
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are  satisfied.   S1(X)  are  given  functions  of  X  and  represent  the  bounds 
on  the  control  u.(t).   The  functions  f (X,u) ,  A(X) ,  and  B(X)  were  pre- 
viously given  in  Section  2.2.   §  is  the  error  in  meeting  the  terminal 
values  of  the  state  variables. 

3.2.   Bounds  on  the  Controls 

The  control  variable  u   is  the  muscle  torque  exerted  at  the 
shoulder  joint  and  u   is  that  exerted  at  the  hip.   For  the  individual 
being  modeled,  the  functions  u   and  u  will  have  upper  and  lower  limits 
which  are  functions  of  the  state  X. 

Samras  [16]  experimentally  determined  the  maximum  muscle  torques 
at  the  shoulder  and  hip  joints  for  various  limb  angles  at  the  joints. 
This  was  done  for  the  same  subject  modeled  in  the  present  study.   These 
measurements  were  made  under  static  conditions  and  the  maximums  in 
either  flexion  or  extension  were  measured  for  the  shoulder  torque  for 
various  values  of  8  and  the  hip  torque  for  various  values  of  i|t.   The  ex- 
perimental bounds  on  the  shoulder  torque  were  then  fitted  by  polynomials 
in  9.   The  hip  torque  bounds  were  expressed  in  polynomials  in  i|i. 

Even  though  each  of  these  bounds  might  be  expected  to  depend  to 

some  degree  on  all  four  state  variables  X3,  X^ ,  X5,  and  XgI  the  bounds 

on  the  shoulder  torque  u  depend  primarily  on  X   and  the  bounds  on  the 

hip  torque  u   depend  primarily  on  X  .   The  measurements  of  Samras  do  not 
2t  0 

include  the  rate  dependence  X^  and  X_.   Although  the  rate  effect  appears 

4      6 

to  be  measurable,  it  is  a  second-order  effect  and  quite  difficult  to 
obtain.   The  control  limit  functions  are  given  in  Figure  8.   These  func- 
tions are  correct  only  for  a  certain  range  of  values  of  the  angles 
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Figure  8.   Unmodified  Control  Limit  Functions 
(Samras  [16]). 


38 

X_  and  X„.   The  values  of  S.  and  S   can  never  be  positive  and  those  of 
3      5  12 

S?  and  s|?  can  never  be  negative.   Whenever  these  sign  conditions  are 

violated  by  extreme  values  of  the  states,  S.  is  set  equal  to  zero.   Also, 

2 

from  extrapolated  measurement  data  an  upper  limit  has  been  set  for  S   at 

160.0  ft-lb  and  a  lower  limit  has  been  set  for  S   at  -100.0  ft-lb. 

3.3.   Torsional  Springs  in  the  Shoulder  and  Hip  Joints 

Our  dynamical  model  and  the  control  limit  functions  of  the 
shoulder  and  the  hip  do  not  account  for  the  stiffness  of  the  shoulder 
and  the  hip  joints  at  the  extremities  of  shoulder  and  leg  movements. 
It  has  been  observed  that  the  shoulder  joints  produce  a  resistance  to 
raising  the  arm  beyond  an  angle  of  8  «  30° .   The  hip  joints  resist  move- 
ment for  i|r  >  120°,  or  for  i|j  <  -35°.   The  effects  of  these  "stops"  are 
important  and  must  be  included  in  the  model,  since  the  film  data  showed 
that  these  limits  were  reached.   There  are  no  data  available  for  the 
stiffness  of  these  joint  stops.   It  was  observed  that,  although  the 
joints  were  not  rigid,  they  were  quite  stiff.   It  was  therefore  decided 
to  use  stiff  torsional  spring  models  at  the  model's  shoulder  and  the 
hip  joints.   These  would  be  active  when  the  stop  angles  were  exceeded. 
For  the  shoulder  the  spring  is  active  for  8  ^  0.5  radian.   For  the  hip 
joint  the  spring  is  active  for  i|i  ^  -0.6  radian  and  i|r  >  2.1  radians. 
The  springs  have  equal  stiffnesses.   One  generates  a  100  ft-lb  torque 
at  the  shoulder  for  a  deflection  of  0.1  radian.   This  corresponds  to  a 
joint  stop  torque  of  the  order  of  the  maximum  voluntary  torque  avail- 
able at  the  shoulder  for  the  deflection  of  0.1  radian.   This  gives  a  spring 

constant  of  K  =  1000  ft-lb/rad.   The  spring  forces  at  the  shoulders 
s 
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would  therefore  be  equal  to  -K  (9-0.5)  for  0  <   0.5  radian  and  those  at 

s 

the  hip  joints  would  be  -K  (\|;-2.1)  for  if   ^  2.1  radians  and  -K  (\|i+0.6) 
for  f  £  -0.6  radian.   These  torques  at  the  shoulder  and  the  hip  joints 
were  added  to  the  voluntary  control  torques  u   and  u  when  the  stops 
were  activated. 


3.4.   Boundary  Conditions 

The  boundary  conditions  for  the  kip-up  maneuver  were  chosen 
from  the  experimental  data  of  Section  2.5.   The  initial  values  selected 
correspond  to  motion  which  has  already  begun.   This  is  beyond  the 
initial  unsymmetrical  motion  which  occurs  on  beginning  the  first  swing. 
This  motion  is  difficult  to  model  and  is  not  important  in  this  basic 
research.   The  final  values  of  the  state  variables  represent  the  model 
atop  the  horizontal  bar  still  moving  upward  and  just  before  body  con- 
tact with  the  bar.   (Once  the  torso  contacts  the  bar,  the  model  is  no 
longer  valid. )   The  actual  motion  in  the  experiment  terminated  shortly 
after  this  point  when  the  gymnast  used  the  impact  of  the  horizontal  bar 
with  his  body  to  stop  himself. 

The  initial  and  final  values  of  the  state  variables  for  the 
optimization  problem  are  listed  in  Table  1. 
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TABLE  1 
BOUNDARY  CONDITIONS  FOR  A  MINIMUM  TIME  KIP-UP  MOTION 

State  Variables       Initial  Value  Final  Value 

CD  X1  =   0.340  X^  =  -2.84 

T  o  f 

Cp  X2  =  -2 .  30  X2  =  -7  .  05 

T  o  I 

3  3 

6  X  =  0.305  Xx  =   2.88 

o  i 

9  X4  =  -0.660  xt  =      0.163 

o  I 

llr  X5  =  -0.087  X^  =   0.436 

T  o  f 

A  X6  =  -1.20  X^  =   0.108 

T  o  f 


3.5.   The  Necessary  Conditions  for  Time  Optimal  Control 

In  this  chapter,  we  look  into  the  necessary  conditions  for  the 
minimum  time  problem  formulated  in  the  previous  chapter.   The  necessary 
conditions  for  optimality  of  motion  for  the  case  when  the  constraints 
on  the  control  are  not  a  function  of  the  states  are  given  in  Reference 
17.   For  the  case  where  control  constraints  depend  on  the  states,  the  nec- 
essary condition  requires  a  modification  in  the  adjoint  equations. 
These  are  obtained  through  a  calculus  of  variations  approach  [18]. 
This  approach  is  used  in  the  following  developments. 
Writing  the  state  equations  of  our  system  as 

X  =  f(X,u)  =  A(X)  +  B(X)u  (3.5.1) 
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we  can  construct  the  cost  function  as 

*f 

J  =  J     1  dt  (3.5.2) 

o 

where  t   is  free.   The  Hamiltonian  is  then  given  by 

H(X,u,X)  =  1  +  \Tf  (3.5.3a) 


=  1  +  \TA(X)  +  \TB(X)u 


(3.5.3b) 


=  1  +  \TA(X)  +  V^tt)^  +  XTB2(X)u2  (3.5.3c) 

where  X(t)  is  the  time-dependent  six-dimensional  column  vector  of 
adjoint  variables.   A,  B,  B.  ,  B_,  and  f  are  the  quantities  as  defined 
by  Equations  (2.2.15)  -  (2.2.22)  in  Section  2.2. 

The  minimum-time  control  policy  u  (t)  will  be  given  by  the  one 
that  minimizes  the  Hamiltonian  (3.5.3),   provided  no  singular  arcs  are 
present.   We  note  that  in  this  case  the  Hamiltonian  is  a  linear  func- 
tion of  the  control  u  and  therefore  the  minimum  with  respect  to  u  occurs 
only  in  the  upper  and  the  lower  bounds  of  u  if  there  is  no  singular 
solution.   Thus,  we  have,  recalling  the  definitions  of  S.  in  Section  3.2, 

T 

(1)  If   X  B  (X)  >  0,  u  =  the  minimum  allowable  value  of  u 

=  S*(X)  (3.5.4a) 

T 

(2)  If   X  B  (X)  <  0,  u  =  the  maximum  allowable  value  of  u 

=  sJ(X)  (3.5,4b) 

T 

(3)  If   X  B  (X)  >  0,  u  =  the  minimum  allowable  value  of  u 

-  -2  -        2  « 

=  S^(X)  (3.5.5a) 

T 

(4)  If   X  B  (X)  <  0,  u  =  the  maximum  allowable  value  of  u 

—  —2  —        2  £ 

=    S^(X)  (3.5.5b) 
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(5)       If      X/V  (X)    = 


X(X)    =    OA 

„(X)   =  oj        x     ' 


=   possible   singular   control.  (3.5.6) 

.   X  y  __v  y,  J-        A 


°r   i  §2  • 

u   and  u  will  be  determined  by  investigating  whether  or  not  there  is 
a  singular  solution  with  respect  to  these  variables. 

The  adjoint  equations  will  be  different  for  the  portions  of  the 
trajectories  for  u   corresponding  to  constrained  and  unconstrained  arcs. 
The  adjoint  equations  are,  in  general,  given  by 

\T  =  -H,  .  (3.5.7) 

This  yields 

(a)  When  neither  u  nor  u  lie  on  a  constraint 

\T  =  -H,  =  -\TA,  -  \\         u   -  XTB    u  .  (3.5.8) 

'X   "  "'X   -  -l,j  1    "  -2>x  2 

(b)  When  any  one  or  both  u   and  u  ,  denoted  by  u.  (i  =  1  or  2),  lie  on 

J.         £  1 

a  constraint  denoted  by  SJ  (j  =  1  or  2)  the  right  side  of  Equation  (3.5.8) 

of  the  adjoint  variables  has  the  additional  term 

T    i 
-  \   B.  SJ 

and  the   equation  can  be  written  as 

2 

\T  =   -    XTA,y  -   XTB  u     -    XTB  u     -      E     5      X  B     S3 

-   ~   X   "  1'X  1      2'X  z        i=l  'X 


(3.5.9) 


where 


6   =  0  if   XTB.  =  0 
i  -  -l 

6.  =  1  if  XTB.  4   0. 
1         -  -1 


The  boundary  conditions  on  the  state  and  the  adjoint  variables 


are 
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X(0)   =  given  =  X         X(0)   =  free 
-  -o        - 

X(tf)  =  given  =  Xf        X(tf)  =  free 


(3.5.10) 


and 


H(t  )  =  (1  +  \Tf)    ^  0  .  (3.5.11) 

1  "  Xf 

The  state  and  adjoint  equations  together  with  the  control  laws 
and  the  boundary  conditions  written  above  form  a  two-point  boundary 
value  problem  (TPBVP)  in  the  state  and  adjoint  variables.   If  these 
equations  can  be  solved,  the  optimal  control,  u  ,  will  be  immediately 
obtained  from  Equations  (3.5.4a)  -  (3.5.6). 

Investigation  of  Singular  Solutions 

T         o 
It  has  been  noted  that  if  \  B  ^  0,  u   cannot  be  determined 

from  the  requirement  that  the  Hamiltonian  is  to  be  minimized  with 

o       T 
respect  to  u  .   The  same  is  true  for  u  when-X  B  =  0. 
X  ^  — '  £ 

Since  the  treatment  for  u   is  the  'same  as  for  u  ,  we  shall 

investigate  a  singular  control  for  only  u  . 

T 
If  the  quantity  \  EL  =  0   only  for  a  single  instant  of  time, 

then  the  situation  is  not  of  much  concern  because  the  duration  of  the 

interval  is  not  finite  and  we  can  simply  choose  u  =  u  (t  )  or  u  (t  ) 

or  0,  where  u  (t~)  =  control  at  the  instant  preceding  t,  u  (t  )  is  the 

instant  exactly  after  t.   The  situation  needs  special  attention  when 

T 
X  B1  =  0   for  a  finite  interval  of  time. 

If   t  ^  t  £  t   is  an  interval  for  which  u   is  singular,  it  is 

clear  that,  for  our  system 

XV  =0     for     t  £  t  £  t  (3.5.12) 
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and  therefore, 


±   (X^)  =  0     for     *!  **  *  *2  (3.5.13) 


or 

T      T 


X  Bx  +  X1^  =  0   for     tt  <  t  <  t2  (3.5.14) 

or,  for  the  interval  t  ^  t  ^  t   the  following  results  must  hold: 

Case  1 

Only  u   is  singular.   u   is  nonsingular.   Since  u   is  not 
1  2  » 

singular,  u   is  on  a  constraint  boundary  and  is  given  by 

u  =  SJ      i  =  1  for  the  lower  constraint 
2    2     J 

j  =  2  for  the  upper  constraint. 
The  adjoint  equations  are  given  by  Equation  (3.5.9) 

^=_XT       T  T      S;J-_XTB2sJ  (3.5.15) 

-  -  x   -  -i»x  1  2-x  2    2    z'x 

'X 

=  B,    [A  +  B,  U.  +  B„  sj]  .  (3.5.16) 

-1     -    -11    -2   2 


From  Equations    (3.5.14),    (3.5.15),    and    (3.5.16),    we  obtain 

[-XTA,       -    XV        S^-XTBSJ2      ]B    +XTB      (A  +    B      S3)    =    0    .  (3. 
'X         ~      2»x      2      -    -2    2,x    -1       -    -1>X             .4      ^ 

It  is  to  be  observed  that  the  necessary  condition  (3.5.17)  is  not 


.1 


explicit  in  u  . 

Case  2 

Both  u   and  u   are  singular.   The  value  of  u   is  no  longer  S 
12  * 

and  the  term   XTB  S^    in  the  X  Equation  (3.5.15)  drops  out  in  this 
case.   Accordingly,  one  obtains 

-XT[A,Y  -  Bn    A]  -  XT[B_    B   -  B     B  ]  u   =  0  .  (3.5.18) 

-  X    -1   -    -  -2,x  -1    -l,x  2   2 
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Proceeding  from  the  assumption  that  u  is  singular,  one  would  also  get, 
for  this  case,  when  both  vl    and  u  are  singular 

-\T[A,„  -  B_   A]  -  XT[B    B  -  B    B  ]u  =  0  .  (3.5.19) 

-   -  X   -2,^  -    -   -1,^  -2    -2,^  -1   1 

From  Equations  (3.5.17),  (3.5.18),  and  (3.5.19)  we  can  see  that 

T        T 
only  if  both  X  B1  and  X  B  are  zero  simultaneously,  is  it  possible  to 

find  a  singular  solution  by  suitable  choices  of  u   and  u   from  the 

T        T         T 
condition  (3.5.13).   If  only  one  of  X  B  and  X  B  ,  say  X  B  ,  is  zero, 

d    T 
the  requirement  —  (X  B  )  =  0  does  not  yield  an  equation  explicit  in 

d     T 
u,  as  observed  in  Equation  (3.5.17).   It  is  thus  required  that  — -  (X  B1 ) 

1  dt   "  _1 

-   0,  which  will  be  explicit  in  u  ,  during  the  interval  t  £  t  £  t 

together  with  the  requirement  that  the  relation  (3.5.17)  is  satisfied 

T 
at  t  =  t  .   These  two  conditions  will  ensure  that  X  B..  =  0  in  the 

interval  t  ^  t  ^  t  . 

It  is  to  be  noted  that  singular  control  computed  by  the  above 

procedure  has  not  been  proved  to  be  the  minimizing  control.   Additional 

necessary  conditions  analogous  to  the  convexity  conditions  for  singular 

controls  have  been  obtained  by  Tait  [19]  and  Kelley,  Kopp  and  Moyer  [20] 

for  scalar  control  and  by  Robbins  [21]  and  Goh  [22]  for  vector  control. 

For  the  general  case  of  vector  control  these  conditions,  summarized  by 

Jacobson  [23],  may  be  stated  as  on  singular  subarcs: 


a 

cTT 


— Ti   H , 

Ldtq   \ 


=  0    if  q  is  odd  (3.5.20) 


and 

,2p 


(-DP  JMAr-  H«  I  *  °  *  (3.5.21) 

™  Ldt2p   UJ 
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d2p 

In  these  equations,  — «~"  H>  (X>X)  is  the   lowest  order  time  derivative 

dt  P   - 

of  H,   in  which  the  control  u  appears  explicitly,  and  q  <  2p. 
u 

For  a  scalar  control,  Equation  (3.5.20)  is  satisfied  iden- 
tically. 

Equations  (3.5.20)  and  (3.5.21)  also  do  not  constitute  suf- 
ficiency conditions  for  minimality.   A  complete  set  of  sufficiency 
conditions  for  singular  arcs  has  not  yet  been  established  in  the 
literature  of  optimal  control  theory  for  a  general  nonlinear  system. 

We  can  see  that  there  are  quite  severe  restrictions  on  the 
existence  of  singular  arcs  in  the  human  motion  problem.   In  the 
numerical  methods  used  in  the  present  work  to  determine  the  optimal 
solution,  only  in  the  method  of  quasilinearization  is  it  necessary  to 
express  the  control  (its  optimal  value)  in  terms  of  the  state  and 
adjoint  variables,  while  in  the  gradient  methods  where  successive 
improvements  are  made  in  the  control  variables,  this  is  not  so.   In 
the  attempts  with  the  quasilinearization  method,  singular  solutions 
were  not  considered  in  the  construction  of  the  two-point  boundary  value 
problem  in  the  state  and  adjoint  variables.   It  was  decided  that  if 
a  solution  to  the  TPBVP  was  obtained  by  quasilinearization,  singular 
arcs  would  be  looked  for  later.   The  gradient  methods  exhibit  singular 
arcs  automatically  if  there  are  any.   The  additional  necessary  condi- 
tions for  singular  arcs  should  be  checked  when  off -constraint  arcs  are 
exhibited  by  the  gradient  method. 
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3.6.  _ The  Solution  Methods 

The  optimal  control  problem  formulated  in  the  preceding  section 
cannot  be  solved  in  closed  form.   Numerical  methods  must  therefore  be 
used  to  find  its  solution.   In  the  optimal  control  theory  literature 
several  numerical  methods  have  been  proposed  for  solving  the  differ- 
ential equations  and  the  optimality  conditions  that  arise  out  of  optimal 
control  problems  such  as  the  present  one.   None  of  these  methods  guaran- 
tees that  a  solution  will  be  obtained  readily,  while  some  of  the  methods 
do  not  guarantee  that  a  solution  may  be  obtained  at  all.   The  methods 
are  all  iterative,  necessitating  the  use  of  high-speed  computers  for 
all  nontrivial  problems.   A  nominal  guessed  trajectory  is  improved 
iteratively  until  the  improved  solution  satisfactorily  meets  all  the 
necessary  conditions. 

Depending  on  whether  the  method  requires  finding  the  first  or 
both  first  and  second  derivatives  of  the  system  equations  with  respect 
to  the  state  and  control  variables,  these  methods  are  called  First-Order 
or  Second-Order  methods,  respectively.   This  is  so  because  they,  in 
effect,  make  first-order  or  second-order  approximations  of  the  system 
equations  with  respect  to  the  state  and  control  variables.   The  first- 
order  methods,  in  general,  have  the  property  that  they  can  start  from 
a  poor  guess  and  make  fast  improvements  in  the  beginning.   They  need 
fewer  computations  in  each  iteration.   But  their  performance  is  not 
good  near  the  optimal  solution  where  the  convergence  rate  becomes  very 
poor.   The  second-order  methods,  on  the  other  hand,  need  a  good  ini- 
tial guess  to  be  able  to  start  but  have  excellent  convergence  prop- 
erties near  the  optimal  solution.   Because  the  second-order  methods 
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need  computation  of  the  second  derivative  of  the  system  equations, 
they  need  more  computing  time  per  iteration,  which  may  be  excessive 
for  some  problems. 

Apart  from  the  first-  and  second-order  methods  mentioned  above, 
there  is  another  class  of  methods  which  tries  to  combine  the  advantages 
of  both  of  these  methods  while  eliminating  the  disadvantages  of  both. 
The  Conjugate  Gradient  Method,  Parallel  Tangent  Method,  and  the  Davidon- 
Fletcher-Powell  Method  fall  into  this  class.   These  methods  work  very 
much  like  the  first-order  method  except  that,  in  the  first-order 
expansion,  the  coefficients  of  the  first-order  term,  or  the  gradient 
term,  is  modified  by  some  transformations.   These  transformations  are 
generated  from  the  modified  gradient  term  of  the  previous  iteration  and 
the  gradient  term  of  the  current  iteration.   This  has  the  effect  of 
using  the  information  that  is  obtained  from  a  second  derivative. 

It  is  not  known  which  of  the  several  methods  used  for  solving 
optimal  control  problems  is  good  for  a  given  problem  and  one  may  have 
to  try  more  than  one  method  in  order  to  obtain  the  solution.   In  the 
published  literature,  most  of  the  illustrations  of  these  methods  are 
simple.   In  these  simple  problems  control  or  state  variable  histories 
do  not  have  wide  oscillations  or  the  system  equations  themselves  are 
not  complicated.   This  makes  it  truly  difficult  for  someone  without 
previous  experience  to  decide  upon  the  merits  of  these  methods. 
There  is  no  preference  list,  and  it  seems  certain  that  there  cannot 
be  one  whereby  a  decision  can  be  made  as  to  which  method  should  be 
tried  first  so  that  a  solution  of  a  given  problem  will  be  obtained 
most  efficiently.   In  this  respect,  deciding  upon  a  computing 
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method  for  a  given  problem  is  still  an  art  and  depends  largely  on  the 
previous  experience  of  the  individual  trying  to  solve  the  problem. 

In  the  attempts  to  solve  the  minimum  time  problem,  the  method 
of  quasilinearization  was  taken  up  first.   This  choice  was  based  on 
several  factors.   This  is  the  only  method  where  the  two-point  boundary 
value  problem  obtained  from  the  necessary  conditions  of   optimality  is 
solved  directly,  and  this  feature  was  found  very  attractive.   As  a  start- 
ing guess,  this  method  requires  the  time  histories  of  the  state  and 
adjoint  variables.   Time  histories  of  the  state  variables  were  available 
from  the  experiments.   (It  was  decided  that  if  the  method  was  success- 
ful for  this  guess,  an  arbitrary  and  less  accurate  initial  guess  would 
be  tried  later.)   When  it  converges,  the  method  has  a  quadratic  conver- 
gence rate.   Also,  in  spite  of  its  being  a  well-known  method  for  solv- 
ing nonlinear  two-point  boundary  value  problems  since  it  was  first 
introduced  by  Bellman  and  Kalaba  [24] ,  its  applications  in  solving 
optimal  control  problems  have  been  very  few.   There  was  thus  an  added 
incentive  for  using  this  method — to  determine  its  usefulness  in  solv- 
ing fairly  complicated  optimal  control  problems. 

Sylvester  and  Meyer  [25]  proposed,  with  demonstrations,  an 
efficient  scheme  for  solving  a  nonlinear  TPBVP  using  the  method  of 
quasilinearization.   This  scheme  was  available  in  the  IBM  SHARE  program 
ABS  QUASI  and  was  used  by  Boykin  and  Sierakowski  [26],  who  reported 
excellent  convergence  properties  of  the  scheme  for  some  structural 
optimization  problems. 

With  this  record  of  success,  the  program  QUASI  was  taken  up  for 
our  problem.   But  with  our  problem  several  difficulties  were  encountered 
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from  the  very  beginning.   First,  the  bang-bang  control  law  obtained 
from  the  necessary  conditions  had  to  be  replaced  by  a  suitably  steep 
saturation  type  control  law.   Second,  a  slight  modification  in  computa- 
tion scheme  was  necessary  when  it  was  found  that  the  method  was  unable 
to  solve  a  simple  example  problem.   The  example  problem  could  be  solved 
with  these  modifications.   But,  in  spite  of  all  these  changes  and  sub- 
sequently, many  attempts  to  generate  a  guess  of  the  adjoint  variables, 
the  method  could  not  be  made  to  work  for  the  human  motion  problem. 
Reasons  for  the  difficulties  encountered  are  discussed  in  detail  in 
Section  3.7. 

During  the  attempts  with  quasilinearization,  it  was  found  that 
computations  of  the  second  derivatives  of  the  system  equations  were 
taking  an   exorbitant  amount  of  time  and  this  was  the  deciding  factor 
for  the  next  choice  of  a  computing  method.   Also,  the  appearance  of  the 
control  function  linearly  in  the  Hamiltonian  put  restrictions  on  the  use 
of  most  of  the  other  second-order  methods. 

The  next  attempts  were  based  on  the  first-order  steepest  descent 
method  proposed  by  Bryson  and  Denham  [27,28].    The  most  attractive 
feature  of  this  method  is  that  the  various  steps  involved  in  it  render 
themselves  to  clear  physical  understanding.   This  method  directly  reduces 
the  cost  function  in  a  systematic  way  and  one  obtains  good  insight  into 
the  basic  steps  in  the  iterative  computations  and  can  make  adjustments 
to  improve  convergence  and/or  stability  with  relative  ease.   These 
features  of  the  method  of  steepest  descent  may  more  than  offset  the 
advantages  of  other  methods  for  some  complicated  problems.   In  the 
attempts  with  this  method,  three  different  formulations  of  the  minimum 
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time  problem  were  tried.   In  the  first  formulation  the  computations  were 
not  pursued  beyond  a  certain  point  due  to  computational  difficulties. 
The  solution  was  obtained  by  the  second  formulation  and  verified  by  the 
third  formulation.   These  attempts  are  discussed  in  Section  3.8. 

3.7.   A  Quasilinearization  Scheme  for  Solving 
the  Minimum -Time  Problem 

In  Section  3.5  the  adjoint  equations  and  the  optimal  control 
laws  (Equations  (3.5.9),  (3. 5. 4a)-(3. 5. 6) )  have  been  derived  for  the 
minimum  time  kip-up  problem.   The  system  equations  and  the  boundary 
conditions  on  the  state  and  the  adjoint  variables  are  given  by  Equa- 
tions (3.5.1)  and  (3.5.10),  respectively.   From  these  equations  we  can 
readily  see  that  if  the  control  variables  U;L  and  u2  appearing  in  the 
system  and  adjoint  equations  are  replaced  by  their  optimal  expressions 
in  terms  of  the  state  and  the  adjoint  variables,  one  obtains  a  non- 
linear TPBVP  in  the  state  and  adjoint  variables.   If  these  equations 
are  solved,  the  optimal  state  and  adjoint  variable  trajectories  will 
be  obtained  and  the  optimal  controls  can  be  constructed  by  using  the 
state  and  adjoint  variables  and  the  optimal  control  laws. 

In  the  TPBVP  in  the  state  and  the  adjoint  variables,  the  final 
time  is  not  a  given  constant  and  is  to  be  determined  from  the  implicit 
relation  (3.5.11).   This  makes  the  problem  one  with  a  variable  end 

point. 

The  method  of  quasilinearization  is  formulated  primarily  for 
a  fixed-end-point  TPBVP.   In  problems  with  variable  end  points,  the 
adjustment  of  the  final  time  is  usually  done  by  a  separate  scheme,  not 
integral  with  the  quasilinearization  scheme.   Long  [29]  proposed  a 


52 

scheme  for  converting  a  variable  end  point  problem  into  a  fixed-end- 
point  problem  with  the  adjustment  of  the  final  time  built  into  the 
quasilinearization  process.   For  the  present  system,  however,  this 
scheme  was  not  practicable  because  the  boundary  condition  (3.5.11) 
becomes  too  complicated  to  handle  in  this  formulation. 

It  was  decided  that  with  a  separate  algorithm  for  adjusting 
the  final  time,  described  later  in  this  section,  the  nonlinear  TPBVP 
with  free  final  time  would  be  converted  to  a  sequence  of  nonlinear 
TPBVP' s  with  fixed  final  times.  Each  of  these  fixed  final  time  problems 
would  then  be  solved  by  the  modified  quasilinearization  algorithm  until 
the  correct  final  time  was  obtained.   The  derivation  of  the  modified 
quasilinearization  algorithm  is  described  below. 

3.7.1  Derivation  of  the  Modified 

Quasilinearization  Algorithm 

The  fixed  final  time  nonlinear  TPBVP  to  be  solved  falls  in  the 

general  class  of  problems  given  by 

dy 

g£=  g(v,t)  .  (3.7.1) 

With  the  boundary  condition 

B«  £(0)  +  Br  £(tf)  +  c  =  0    tf  =  given  (3.7.2) 

V,  g,  and  c  are  of  dimension  n,  B.  and  B  are  matrices  of  dimension 
(n  xn).   It  is  being  assumed  that  the  TPBVP  has  been  defined  for  the 
interval  0  £  t  <   t   for  some  given  t  >  0. 

In  the  state  and  adjoint  equations,  if  the  expressions  for  optimal 
control  in  terms  of  the  state  and  the  adjoint  variables  are  used  for  the 
control  variables,  one  obtains 
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X  =  f(X,u°(X,X))  =  F(X,X)      (say) 


(3.7.3a) 


and 


X  =  -h£(X,u°(X,X),X)  =  G(X,X)   (say).  (3.7.3b) 

In  the  formulations  of  the  TPBVP  given  by  Equations  (3.7.1)  and 


(3.7.2),  it  may  be  seen  that  for  the  kip-up  system,  n=12, 


and 


X 

_X  _ 

I  = 

F 
G 

"$f  " 

X 
-o 

• 

0  0 

1  o 


B  = 


I   0 
0   0 


c  = 


The  0  and  I  appearing  in  the  matrices  B^  and  Br  represent  6x6  order  null 
and  unit  matrices,  respectively. 

Let  z(t)  be  an  initial  guess  vector  for  y(t)  which  satisfies 
the  boundary  conditions  (3.7.2).   If  g(y,t)  is  approximated  by  its 
Taylor  series  expansion  about  g(z,t),  keeping  only  the  first-order  term, 
one  obtains 


g(y,t)  =  g(z,t)  + 


Sg 


oy 


(y-z) 


y=z 


Let 


W  = 


w.  . 
1J 


-  y=z 


dgi 


,    so  that, 


(3.7.4) 


y=z 


th 


or  W     partial  derivative  of  the  i   element  of  g  with  respect  to  the 
ij 

j   element  of  y,  evaluated  at  y  =  z. 


or, 


where, 
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With  the  above  approximation  of  g(y>t),  Equation  (3.7.2)  becomes 

dy 

—  =  g(x,t)  +  W(z,t)(y-z) 
at   -  -       -    -  - 


de     dz 

"=-_!+  g(z,t)  +  W(z,t)   e 
at     at   -  -       -     - 


e  =  y(t)  -  z(t)  =  error  in  the  guess  z(t). 

Rearranging  the  above  equation,  one  obtains 

de  dz 

-1  -  W(z,t)€  =  -  ^1+  g(z,t).  (3.7.5) 

Since  z(t)  is  chosen  to  satisfy  the  boundary  conditions, 

B.   z(0)  +  B  z(t_)  +  c  =  0. 

JC   —        r  —   I     -    — 

Subtracting  this  equation  from  Equation  (3.7.2),  one  obtains  the 

boundary  conditions  on  the  error  e(t)  as 

B.  e(0)  +  B   e(t„)  =  0.  (3.7.6) 

Hj  —  r  —   1     — 

Equations  (3.7.5)  and  (3.7.6)  form  a  linear  TPBVP  in  e(t)  which,  when 

solved,  will  give  the  values  of  the  error  between  the  guessed  solution 

z(t)  and  the  actual  solution  y(t)  based  on  the  linearized  expressions 

of  the  right  side  of  Equations  (3.7.1)  about  the  guessed  solution  z(t). 

Because  of  using  the  linearized  equation  instead  of  the  full  nonlinear 

equations,  the  values  of  e(t)  obtained  by  solving  Equations  (3.7.5)  and 

(3.7.6)  will  not  be  the  actual  error  between  the  guess  z(t)  and  the 

solution.   However,  a  new  guess  of  y(t)  will  be  obtained  from  e(t)  by 

z'(t)  =  z(t)  +  7]  e(t),     0  <  T]  <  1.  (3.7.7) 
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The  algorithm  of  Sylvester  and  Meyer  uses  11  =  1  for  all  the 
time,  which  is  the  usual  quasilinearization  algorithm.   It  was  found, 
while  solving  a  simple  example  problem,  that  without  the  incorporation 
of  a  multiplier  7]  in  the  expression  (3.7.7),  i.e.,  using 

z'(t)  =  z(t)  +  e(t) 

the  method  was  unstable.   The  convergence  property  of  the  scheme  with 
the  incorporation  of  the  multiplier  7]  can  be  understood  for  small  values 
of  T)  by  comparison  with  the  step-size  adjustment  procedure  of  the  usual 
steepest  descent  algorithms.   The  mathematical  proof  for  the  convergence 
property  follows  the  proof  of  Miele  and  Iyer  [30]  and  is  now  given. 

The  integral  squared  norm  of  the  error  in  the  guessed  solution 
z(t)  can  be  expressed  by  the  integral 
t. 


J  =  f  '  {z  -  g(z,t)}T  [z  -  g(z,t)}  dt. 
Similarly,  the  error  in  the  solution  z'  (t)  =  z(t)  +  7]  e(t)  is  given  by 


** 


J 


'  =  !    i*   -  !<*'.*>}    ik'  -  s(z',t)}  dt 


If  71  is  sufficiently  small,  one  can  write 
g(z',t)  =  g(z,t)  +  g,z(z,t)  71  e_ 


where   g,   =  g.  I    =  W(z,t)  (from  Equation  (3.7.4)). 
*-  y=z 


Also,  for  all  values  of  7], 
z  =  z  +  7]  e  . 
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From  these  results,  one  obtains,  for  small  values  of  T) , 

j'  -  j  =  21)  /   (4  -  g(z,t)}T  [e  -  we}  . 
o 

Since  e(t)  satisfies  the  differential  equation  (3.7.5),  this  finally 

yields: 

*f 
J7  -  J  =  -  2T}  J   ||  z  -  g(z,t)J|  dt 
o 

=  a  negative  quantity. 
Thus,  for  sufficiently  small  values  of  T),  the  reduction  in  the  cost  is 
assured. 

In  the  quasilinearization  algorithm  z  (t)  takes  the  role  of  z(t) 
as  the  new  guess  of  y(t)  and  the  process  is  continued  until  the  error  in 
satisfying  the  differential  equations  is  reduced  to  an  acceptable  value. 

The  linear  TPBVP  of  the  error  Equations  (3.7.5)  and  (3.7.6)  is 
solved  as  follows: 

The  time  interval  t  =  0  to  t  =  t   is  divided  into  m  small  inter- 
vals.  This  results  in  m+  1  values  of  t  at  which  the  solution  will  be 
computed.   Equation  (3.7.5)  can  be  written  in  a  finite  central  differ- 
ence scheme  as 

°i+i - gj  „ A +  2-i+i   VW]  -ei +  Jul    Si+i-Si  e r-ui  +  z-i  Vi+Vi 

—  h < 2 >— J   2 = h" +iV— 2 '  —2 ) 

i  i 

th 
where  h  =  t    -  t  .   The  subscript  i  denotes  values  at  the  i 
i    i+1    i 

station,  i = 1,2, . . . ,m+l.   Rearranging  and  simplifying  the  above  expres- 
sion one  obtains 

(I  +  i  h.W.)e.  +  (-I  +  i  h.W.)e.   =  r.  (3.7.8) 

2   i  i  -i       2   l  l  -l+l    -i 

i  -   1,2, ...  ,m 
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where 


(  -1+1   -1  1+1    u 

r.  =  z .   -  z  .  -  h .  g  I , ) 

-l   -l+l  -l   i  -  V   2  2   y 

/Z.  .,  +  Z  .     t.  .,+t.N 


(3.7.9) 


(3.7.10) 


and 

I  =  unit  matrix  of  dimension  n  xn  . 
The  boundary  conditions,  Equation  (3.7.2),  reduce  to 


B.  e  +  B  e    =  0 
Z  -1  r  -m+1 


(3.7.11) 


Equation  (3.7.8)  can  be  cast  into  the  following  convenient  recursive 


expression 


where 


and 


e.  +  D.  e.  „  =  s. 
-l    i  -i+I   -l 


D.  =  (I  +  Ih.  W.)-   (-1+  \       h.  W.) 


*   1        -1 
s.  =  (I  +  -  h.  W  )    r 
-l       2   l   l     -i  J 


(3.7.12) 


(3.7.13) 


By  repeated  substitution,  equation  (3.7.12)  yields  the  following  rela- 
tionship between  e.,  and  e 

-1     -m+1 


m 


m 


where 


£l  =  T+  (-1)   (TTD.)  z 
i=l 


m     ._1   i-1 
T  =  s  +   S  (-1)1    (  TT  D  )s 
i=2         j=l  J 


(3.7.14) 


(3.7.15) 


or,  multiplying  by  B«  on  both  sides  of  Equation  (3.7.14), 


.  "1 


m 


Bih  "  Bzl +  (-1}    h  J^i  -Vi 


(3.7.16) 


Equations  (3.6.16)  and  (3.7.11)  can  be  solved  simultaneously  for  e 


and  G   ,  to  give 

-nvt-1 
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6   ,  =  -C"1  B,  T  (3.7.17) 

-m+1        Xi  - 

where 

m 

C  =  (-I)1"  B.  (  It   D.)  +  B   .  (3.7.18) 

x  .  ^  i     r 
i=l 

With  e  „  determined,  s  ,  e  H . . . . , e„  are  determined  in  succession  by 
-m+1  -m  -m-1     -1 

using  the  recursive  relation  (3.7.12).   Then,  the  new  iterate  is  given 
by 

zi+1(t)  =  zX(t)  +  7]  e(t). 

The  stopping  condition  of  the  algorithm  is  given  by  the  fact  that 

e   should  be  small.   When  thev  are  small,  it  may  be  seen  from  Equation 
~x 

(3.7.8)  that  the  quantities  r.  will  also  be  small.   From  Equation  (3.7.9) 

it  is  also  seen  that  small  r.  correspond  to  satisfying  the  central  differ- 

-l 

ence  expression  of  the  differential  equations.   That  is,  the  finite  dif- 
ference equation  error  must  be  small.   This  does  not  mean  that  the  dif- 
ferential equation  error  is  small  unless  the  intervals  for  the  differences 
are  sufficiently  small. 

The  IBM  SHARE  program  ABS  QUASI  is  a  program  of  the  procedure 
outlined  above  without  the  provision  of  the  multiplier  T|  in  Equation 
(3.7.7),  and  therefore  is  for  T|  =  1.   The  program  was  modified  to  intro- 
duce and  to  adjust  T]  to  get  the  desired  convergence.   The  algorithm  may 
be  described  by  the  following  steps: 

1.   Set  up  the  matrices  B,  and  B  and  guess  a  nominal  trajectory 

X,  r 

z(t)  that  satisfies  the  boundary  conditions  (3.7.2). 
Set   ITER  =  0     (ITER  for  Iteration) 
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2.  Do  the  following  for  i  =  l,2,...,m. 

a   Find  r  and  W  as  defined  in  Equations  (3.7.9)  and 
-i      i 

(3.7.10).   Find  the  largest  element  of  r.,  searching 
between  the  elements  of  each  r.,  for  i =  1,2, . . . ,m. 
Call  it  E2MAX.   If  E2MAX  <  specified  maximum  allowable 
error,  print  out  z.  and  stop  the  computations, 
b.   Using  Equation  (3.7.13),  find  D  and  s±. 

3.  Calculate  T  and  C  according  to  Equations  (3.7.15)  and  (3.7.18). 
Calculate  the  integral  norm  of  the  error 

i=2 
(here  the  norm  is  defined  by  the  sum  of  squares  of  the  elements 

of  the  vector  r.).   Set  Jl  =  J2. 

-l 

4   Find  e  .,  using  Equation  (3.7.17). 
-m+1 

5.  Decide  upon  a  value  of  T|.   Discussions  on  the  choice  of  11  will 
be  presented  in  a  later  section. 

6.  Generate  and  store  e  ,  e    ,  e    ,...,e  using  Equation  (3.7.12) 

-m  -m-i   -m-z     —x 

Generate  the  new  guess  z.,  i =  1 ,2, . . . ,m+l,  by  doing 

z.  =  z.  +  T)  e. . 
-l   -l     -i 

Set   ITER  =  ITER+  1. 

7.  Find  J2  and  r  ,  i  =  l,2,...,m  and  find  E2MAX  as  in  step  2. 

8.  If  J2  >  Jl  (unstable),  go  to  step  11. 

9.  Stop   if  E2MAX  <  a  prescribed  value. 

10.  If  J2  <  Jl  to  to  step  12  to  continue  to  the  next  iteration. 

11.  If  this  step  has  been  performed  more  than  a  specified  number 
of  times  in  this  situation,  go  to  step  13.   If  not,  store  the 
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value  of  the  current  T)  and  J2.   Recover  the  values  of 

z  of  the  previous  iteration  by  doing 
-i 

z  «•  =  z .  -7)  e .  . 
-l   -l     -l 

Reduce  the  value  of  T).   Generate  the  new  a  by  doing 

z.  =  z.  +  7]  e.  .  V 

-i   -i     -i 

y 
Go  to  step  7. 

12.  If  this  step  has  been  performed  more  than  a  specified  number 
of  times  or  if  step  11  has  been  performed  at  least  once  in 
this  iteration,  go  to  step  13.   If  not,  store  the  value  of  J2 
and  7).   Recover  the  new  z.  as  in  step  11,  this  time  by  increas- 
ing 7]  to  increase  speed  of  convergence.   Go  to  step  7. 

13.  Find  out  the  minimum  values  of  J2  that  have  been  obtained  in 
steps  11  and  12.   If  this  value  is  greater  than  Jl ,  stop  com- 
putations and  look  for  the  cause  of  the  instability.   If  this 
value  is  less  than  Jl,  recall  the  7]  corresponding  to  this  J2  and 
regenerate  the  z.  for  this  case.   Go  to  step  2. 

In  this  program  the  value  of  7]  was  selected  from  the  consider- 
ation that  at  the  station  m+1  the  maximum  value  of 

should  be  less  than  a  prescribed  value  (j  in  parentheses  represents  the 
jth  element  of  the  vector).   This  procedure  limits  the  changes  in  the 
existing  trajectory  by  limiting  the  magnitude  of  the  maximum  fractional 
change  in  the  terminal  values  of  the  variables  not  specified  at  the 
final  time. 
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Whenever  an  iteration  was  found  unstable,  T)  was  reduced  by 
half.   When  there  was  an  improvement,  a  linear  extrapolation  formula 
was  used  to  increase  the  value  of  T]  so  that  the  norm  of  the  error  J2 
would  decrease  to  a  desired  value.   In  such  an  attempt,  however,  T\ 
was  not  allowed  to  increase  beyond  a  certain  multiple  of  its  existing 
value. 

3.7.2.   Approximations  of  the  Optimal  Control 
for  the  Kip-Up  Problem 

The  method  of  quasilinearization  disregards  the  question  of 

singular  solutions  in  the  present  investigation.   It  was  found  in 

Section  3.5  that  there  were  many  requirements  for  the  existence  of 

a  singular  control  arq  and  the  necessary  conditions  for  the  existence 

of  singular  controls  are  quite  complicated.   It  was  decided  that,  before 

going  into  those  cases,  extrema  without  singular  arc  would  be  looked  for 

first  and  if  the  computational  method  was  successful,  singular  subarcs 

would  be  looked  for  later. 

If  the  singular  solutions  are  not  considered,  the  optimal 

o      o 
controls  u   and  u  become  bang-bang  and  are  given  by  Equations  (3.5.4a) 

J-  £i 

through  (3.5.5b). 

It  is  convenient  to  rewrite  these  equations  at  this  point 
in  the  following  way: 


if  -xV^x) 


<  o, 

o 
Ul   = 

s*qp 

=  o, 

o 
Ul   = 

0 

>  0, 

o 

ul  = 

S^(X) 

(3.7.19) 
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and 


where 


and 


If   -X  B  (X) 


<0,  u°  =  S*(X) 


=  0,   u  =  0 


(3.7.20) 


>  0,  ul   = 


s*<x) 


Let  these  control  laws  be  expressed  by  the  following  expressions: 
u°  =  S1  sgn  (-X^)  (3.7.21) 


i  =  1  or  2 


u°  =  S2  sgn  ("XTB2) 


sgn  (x)  -  sign  of  x 


=  S1  when  (-\TB.)  <  0 
l         -  -l 


=    S2  when  (-\TB.)  >  0 
l         -  -l 


(3.7.22) 


(3.7.23) 


(3.7.24) 


With  these  expressions  for  the  control  laws,  the  state  equations  (3.5.1) 
become 

X  =  A(X)  +  B  •  S   sgn  (-X^)  +  52  •  S2  sgn  (-X  B^  (3.7.25) 


=  F(X,\) 


(cf.  Equation  (3.7.3a)). 


In  the  present  quasilinearization  algorithm  the  derivatives  of 
F(X,X)  with  respect  to  both  X  and  X  in  constructing  the  matrix  W^^ 
(Equation  (3.7.4))  are  needed.   One  can  see  that  computation  of  the 
derivative  of  F(X,X)  with  respect  to  X  will  occur  only  as  a  general- 
ized function  with  no  numerical  value  for  computation  in  the  limit, 
because  X  appears  only  in  the  argument  of  a  sign  function.   Instead 
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of  proceeding  to  the  limit  the  following  approximation  of  the  bang-bang 
control  was  made.   For  i  =  1  and  2 

r  1  T      1 

S  (X)   if   AA. (-\  B.)  <  S. (X) 
l  -         l   -  -l      l  - 

u°  «*  u*  :    AA.(-\TB.)   if   S2(X)  >  AA.(-\TB.)  ^  S1(X)  (3.7.26) 

1--1         l-       1--1      1  - 

S2(X)   if  AA.(-\TB.)  >  S2(X) 
L  1  -  1--1      1  - 

where  AA,  and  AA  are  two  positive  constants.   This  function  u.  of 

T  12 

AA. (-\  B.)  is  called  the  saturation  function  (sat)  when  S.  and  S.  are 

unity.   The  change  made  in  the  optimal  control  is  shown  graphically  in 

o      * 
Figure  9  for  u.  and  u. ,  i  =  1  or  2. 

The  controls  u  and  u.   have  been  plotted  against  the  function 
l      l 

T 
-X  B  near  a  switch  point  in  Figures  9a  and  9b,  respectively.   It  can 

-  -x 

*  o 

be  seen  that  the  approximation  u.  differs  from  the  optimal  control  u. 


ly  in  the  portion  KL.   By  increasing  the  value  of  AA. ,  u.  can  be  made 


on 


o 
closer  to  u. . 

l 


With  the  above  approximation  of  the  bang-bang  control  by  a 

"saturation"  control  represented  by  Equation  (3.7.26),  one  would 

be  able  to  find  the  derivative  of  the  control  u.  (and  hence,  of  F)  with 

respect  to  X.      The  derivative  will  be  zero  when  the  control  is  on  a 

constraint  S  (X)  and  will  be  nonzero  on  the  arc  KL  which  appears  near 
i  - 

a  switching  time.   In  order  to  represent  the  saturation  control  (3.7.26) 
by  its  linearization,  it  is  necessary  that  at  least  one  of  the  W. , 
i  =  l,2,...,m,  which  contains  the  information  of  the  first  derivatives, 
be  computed  on  the  arc  KL.   Otherwise,  this  vital  portion  of  the  control 
would  go  unaccounted  for  in  the  linearized  equations.   This  may  happen 
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if  the  arc  KL  is  too  steep  for  a  given  selection  of  integration  stations. 
In  such  a  case,  when  none  of  the  W.  is  computed  on  the  switching  por- 
tions of  the  control  variables   like  the  arc  KL, the  TPBVP,  Equations 
(3.7.5)  and  (3.7. 6),  cannot  be  solved  as  explained  below. 

First,  from  a  physical  reasoning,  it  can  be  seen  that  the  linear- 
ized state  equations  would  get  decoupled  from  the  adjoint  equations  if 
all  of  the  W. ,  i  =  l,2,...,m,  show  zero  derivatives  of  F  with  respect 
to  \.      This  means  that  the  first  six  equations  of  (3.7.5)  would  get 
decoupled  from  the  last  six.   But  the  boundary  conditions  (3.7.6)  is  such 
that  they  are  for  the  first  six  variables  of  e(t)  only.   Therefore,  this 
results  in  a  situation  where  there  are  six  first-order  equations  with 
twelve  boundary  conditions,  a  situation  which,  in  general,  does  not 

have  a  solution. 

From  the  point  of  view  of  computations  with  the  present  algorithm 
it  can  be  shown  that  the  matrix  C  in  equation  (3.7.17)  cannot  be  used  to 

solve  for  e  ., . 
-m+1 

With  the  system  equations  defined  as  (Equations  (3.7.3a)  and 

(3.7.3b)) 

X  =  F(X,\) 
\  =   G(X,\) 


one  obtains 


[W]. 


2»x    F-'x 


S'x    S'x 


th 


where  the  subscript  i  represents  the  i    interval 
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Let  P  . ,  Q . ,  E . ,  DD . ,  R .  and  M  ,  j 
3         3         3  3         3  3 


1,2,3,  and  4,  represent 


6x6  matrices.   Let 


[W]±  = 


P  P  ~ 

1  2 


P  P 

L    3        *4j. 


where     P_  = 


N 


[0]      for   all   points   other   than  those  lying  on  the 


arcs   such   as  KL   in  Figure  9b.      Then 


1  ,-1 

[I   +    -  h.    W.] 


1   +    \  hi    Pl  \  \   P2 


-1 


1   ,       „ 
2hiP3 


1   +    2   hi   P4 


\ 


Q,- 


°-3        Q4 


It  may   be   seen  that   Q2  =    [0]      if      ?2  =    [0] . 
From  Equation    (3.7.13), 

°i  =  Hh^_1MhiW3 

ft         Q2 


LQ3         Q4 


^1        R2 


LR3        R4J 


El        E2 


LE3        E4J 


(say)  , 


where 


Rl        R2 


R3        R4 


p   +    \  hiWi]     • 


and 


*2   =   I  hi    P2  =    r°]       "      P2  =    [°] 
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It  may  again  be  seen  that 

E2  =  Q1R2  +  Q2R4  =  [°]    SlnCe   R2  =  Q2  =  [°]    Wh6n    P2  =  [°L 

Now,  consider  the  product  D.  •  D    •  D    of  any  three  succes- 

m  . 

sive  D  in  the  product  TT  D..   Let  this  multiplication  be  expressed  as 
i  i=l   1 


Ell    E21 


E31    E41 


E12    E22 


E32    E42 


E13    E23 


E33    E43 


DD_ 


DD. 


DD3     DD4 


The  expression  for  DD  is 

DD2  =  Ell  E12  E23  +  E21  E32  E23  +  Ell  E22  E43  +  E21  E42  E43' 

If  the  upper  right  elements   Egl  =  *22  =   %23   =  [0],  it  may  clearly  be 

seen  that  DD2  =  [0] . 

In  the  product 


m 

l\ 

M2 

TT     D.  = 

M3 

M4 

one  would  therefore  obtain 


M  =  [0]   if   P  =  [0]   for  all  i,   i  =  l,2,...,m. 
2  2 

With   M  ==  [0] ,  the  expression  for  C  (Equation  (3.7.18)) 


m 


m 


D  =  (-1)   B.  TT  D  +  B 
1   i=l   X 


becomes,  after  using  the  values  of  B^  and  B^, 


Go 


BZ  = 


0  0 

1  0 


C  =  (-1) 


m 


I    0 


m,_  m2 


I   0 


0   0 


=    (-1) 


m 


I    0 


M    0 


Ml 


Thus,  C  becomes  singular  if  F,^  =  0   for  all  i,  i  =  l,2,...,m. 

For  a  given  step  size  or  subdivision  in  the  t-interval  (0,t  ) 
the  singularity  of  C  sets  an  upper  limit  on  the  steepness  of  the  arc  KL 
of  Figure  9b  which  can  be  used  for  approximating  the  bang -bang  control. 
This  steepness  of  the  arc  KL  depends  on  7^  (or  AA2)  for  a  given  X(t) 
and  \(t).   So,  in  selecting  AA  and  AA  ,  we  have  to  make  sure  that  they 
are  not  too  large.   When  the  solution  is  obtained,  however,  the  slope  of 
KL  does  not  depend  on  the  choice  of  AA  and  AA2  if  we  decide  to  choose 
AA  =  AA  -   AA,  say.   This  is  because  AA  will  be  "absorbed"  within  X(t) , 
acting  as  a  scale  factor  on  X.   In  fact,  if  we  define  X  as  AA«\,  the  state 
and  adjoint  equations  of  our  system  do  not  change,  and  we  could  therefore 
select  AA  =  AA  =1.   The  slope  of  the  arc  KL  in  the  solution  depends 
on  the  value  of  AA» \  and  not  on  AA  alone.   However,  when  the  solution 
has  not  yet  been  obtained,  the  best  value  of  AA  need  not  be  1. 

If  the  final  time,  t  ,  selected  for  our  problem  is  much  larger 
than  the  minimum  time,  for  a  bang-bang  optimal  solution,  it  can  be 
expected  that  the  arc  KL  will  have  a  relatively  small  slope.   If  the 
final  time  is  gradually  reduced,  this  slope  will  increase  until  it 
becomes  so  large  that  the  matrix  C  becomes  singular  as  explained  above. 
At  that  point  the  solution  obtained  for  the  smallest  tf  will  very 
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nearly  be  a  bang-bang  control  and  will  represent  the  approximate  solu- 
tion of  the  minimum  time  problem. 

The  computation  of  the  optimal  final  time  via  the  quasilinear- 
ization  method  was  done  according  to  the  stopping  condition  outlined 
above.   A  final  time  would  be  guessed,  and  the  TPBVP  (Equations  (3.7.1) 
and  (3.7.2))  would  be  solved.   The  final  time  would  then  be  reduced  by- 
reducing  the  integration  step  size,  and  the  TPBVP  would  again  be  solved. 
The  process  would  be  continued  until  the  matrix  C  becomes  singular  and 
the  TPBVP  could  not  be  solved  any  further. 

3.7.3.   A  Simple  Example  Problem  for 

the  Method  of  Quasilinearization 

A  simple  problem  was  first  taken  up  to  explore  the  various 

features  of  the  quasilinearization  algorithm  constructed  above.   The 

system  was  defined  as 


X1=X2 


X2=U 


(3.7.27) 


X1(0)  =  0 

x2(o)  =  0 

w  - 1 

VV  ■  ° 


(3.7.28) 


The  constraints  on  u  were 
-1  £  u  ^   1 

The  cost  function  to  be  minimized  is  the  final  time  t^. 


(3.7.29) 


The  adjoint  equations  for  the  system  are 
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(3.7.30) 


The  optimal  control  u   is 


u  =  sgn  (~X2) 


(3.7.31) 


With  the  approximation  of  the  optimal  control 


u 


sat  (-AA»\  ) 


(3.7.32) 


the  state  and  adjoint  equations  become 


X1  =  X2 

X  =  sat  (-AA-X  ) 
2  z  (3.7.33) 

\=° 

h   =  "Xl  • 

The  boundary  conditions  are  given  by  Equations  (3.7.28).   The  analytical 
solution  of  the  optimal  control  is  given  by 

-1 


>» 


*f  =  2 

x1  =  -i 

\2=  t-1 


u 


X  =  -1  +  2t  - 1  t  2  \  f  or  1<  t  ^  2 

X  =  2  -  t 
2  J 


~\ 


u  =  1 


X. 


1  .2 


t  \     for   0  £  t  <  1 


2       J 
The  solution  is  shown  graphically  in  Figure  10. 

The  problem  was  solved  numerically  by  solving  Equations  (3.7.33) 
and  (3.7.28)  for  different  final  times  t   by  quasilinearization. 
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Figure  10.   Graphs  of  Optimal  and  Nearly  Optimal 

Solutions  Obtained  via  Quasilinearization 
for  Simple  Example. 
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The  solutions  for  t  =  2.25,  2.05,  and  2.005  were  obtained  and  are  shown 
in  Figures  10a,  b,  and  c,  respectively.   With  t   =  2.00  the  matrix  C 
became  singular  in  iteration  10.   The  theoretical  solution  (t  =2.00) 
is  shown  in  Figure  lOd.   The  results  for  t  =2.005  is  a  reasonably  good 
approximation  of  the  optimal  solution. 

For  these  problems,  25  time  subdivision  intervals  were  used  . 
The  initial  guess  was  deliberately  poor,  as  shown  in  Figure  lOe.   The 
solutions  were  obtained  in  9  to  11  iterations  from  this  guess  in  all 
these  cases. 

The  program  originally  written,  according  to  Sylvester  and 
Meyer  [25]  ,  did  not  have  the  provision  for  the  amount  of  adjustment  T) 
in  the  iterations.   Their  method  was  found  to  be  unstable  for  some  of 
the  initial  guesses. 

3.7.4.   The  Results  With  the  Kip-Up  Problem 

The  kip-up  problem  was  taken  up  after  the  method  of  quasilinear- 
ization  was  found  successful  in  the  case  of  the  example  problem.   In  the 
kip-up  problem,  however,  many  difficulties  were  faced  from  the  very 
beginning  and  the  problem  could  not  finally  be  solved  by  this  method. 
A  major  difference  between  the  human  motion  problem  and  the  example 
problem  or  the  problem  solved  by  Boykin  and  Sierakowski  [26]  is  very 
prominent.   In  the  latter  problems,  the  control  variables  switched  only 
once  from  one  boundary  to  another  in  the  entire  trajectory,  whereas, 
in  the  human  motion  problem,  there  were  many  such  switchings.   This 
made  the  human  motion  problem  less  amenable  to  iterative  methods. 
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The  program  for  the  human  motion  problem  was  extremely  lengthy, 
taking  more  than  eleven  hundred  statements  and  requiring  the  use  of  the 
large  core  of  the  computer  (IBM,  360-65).   The  subroutine  (DIFEQ)  which 
generated  the  right  side  of  the  state  and  the  adjoint  equations  and 
their  derivatives,  turned  out  to  be  quite  lengthy  and  required  an 
exorbitant  amount  of  computing  time.   The  quasilinearization  program 
called  this  subroutine  at  every  station  of  the  total  interval  and  at 
every  iteration.   As  a  result,  the  program  required  a  tremendous  amount 
of  computing  time — about  40  seconds  per  iteration.   All  computations 
were  done  in  double  precision. 

Several  sources  of  difficulty  were  detected  in  the  unsuccessful 
attempts  to  solve  the  human  motion  problem  by  quasilinearization.   The 
central  issue  was  the  tremendous  amount  of  computations  with  accuracy 
required. 

The  total  time  interval  was  divided  into  100  equal  parts  (seg- 
ments) .   As  the  first  (and  the  only  one)  guess  of  the  final  time,  these 
segments  were  chosen  to  be  0.0216  second   each.   This  made  the  final, 
time  (2.16  seconds)  equal  to  the  time  taken  by  the  gymnast  to  do  the 
actual  maneuver.   With  100  segments,  the  program  required  the  large 
core  of  the  computer.   Even  though  a  larger  number  of  smaller  segments 
would  have  been  preferred,  this  could  not  be  done  because  it  was 
intended  not  to  go  beyond  the  large  core.   It  was  later  found  that 
numerical  integration  of  the  nonlinear  state  equations  needed  time 
steps  at  least  as  small  as  one-sixth  of  what  was  taken  for  quasilinear- 
ization.  It  was  initially  hoped  that,  since  the  quasilinearization 
algorithm  solves  the  linearized  equations  instead  of  the  nonlinear 
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equations,  the  central  difference  solution  would  be  stable  for  larger 
integration  step  sizes.   This  was  not  the  case. 

It  was  proved  theoretically  that  with  perfect  precision  the 
method  was  convergent  for  all  initial  guesses.   However,  it  is  well 
known  for  the  unmodified  method  of  quasilinearization  (T|  =  1)  that  the 
method  is  convergent  only  for  certain  initial  guesses.   It  may  therefore 
be  fair  to  expect  that  even  with  the  modified  method,  the  rate  of  con- 
vergence should  depend  on  the  initial  guess  and  for  some  initial  guesses, 
this  convergence  may  be  extremely  poor.   For  this  reason,  several  initial 
guesses  were  tried.   The  guess  of  the  state  variables  was  taken  from  the 
experimental  data  which  were  shown  to  agree  .veil  with  the  computed  motion 
in  Chapter  2.   Different  initial  guesses  for  the  adjoint  variables  were 
tried.   The  first  attempts  simply  used  constant  values  for  all  the 
adjoint  variables.   Convergence  from  this  guess  was  nonexistent.   Next, 
attempts  were  made  to  generate  the  adjoint  variables  from  forward  inte- 
gration of  the  adjoint  equation  with  the  guessed  values  of  the  state 
variables  and  a  guessed  value  of  adjoint  variables  at  time  t  =  0.   In 
these  cases  the  integrations  were  unstable  with  large  numbers  generated. 
The  method  was  not  pursued  further. 

Lastly,  the  method  suggested  by  Miele,  Iyer  and  Well  [31]  was 
tried  for  generating  the  initial  guess  of  the  adjoint  variables.   In 
this  method,  an  auxiliary  optimization  problem  is  constructed  from 
the  original  problem.   It  tries  to  make  an  optimal  choice  of  the  adjoint 
variables  such  that  the  cumulative  norm  of  the  error  in  satisfying  the 
adjoint  equations  for  a  given  state  variable  trajectory  is  minimized. 
This  is  performed  as  follows. 
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Suppose  u  (X,\)  is  the  optimal  control.   The  state  and  adjoint 
equations  may  then  be  .vritten  as 

X  =  f(X,U*(X,>.))  =  F(X,X)  (3.7.34) 

X  =  -&,U*(X,»)X  =  G(X,X)  •  (3.7.35) 

Suppose  we  have  a  guess  of  the  state  and  the  adjoint  variables 
given  by  X(t)  and  X(t) .   Since  X(t)  and  X(t)  do  not  satisfy  Equations 
(3.7.34)  and  (3.7.35) ,  let 

e  =  X  +  Ax,u*(x,X))X 

so  that  e  is  the  error  in  the  adjoint  equations  (3.7.35).   The  above 

equation  can  be  rewritten  as 

T 
X=  e  -  f,x(X,u*(XfX))X  (3.7.36) 

Now,  consider  the  optimal  control  problem,  where 

a.  X(t)  is  a  given  function  of  time  t 

b.  X(t)  is  the  state  variable 

c.  e  is  the  control  variable 

d.  The  system  equation  is  given  by  Equation  (3.7.36) 

e.  The  cost  function  to  be  minimized  is 


1  P  T  ^ 

j  =  —  i  e  e  dt 

2    2  "  -  - 
o 

f.  t   is  the  final  time  (fixed)  of  the  original  problem 

g.  Boundary  conditions  are 

A.(0)   ar.d   >.(t  )  are  free. 
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For  the  optimal  control  problem  posed  above,  we  can  construct 

T 

the  Hamiltonian,  using  the  Lagrange  multipliers  e^  (6-dimensional 

vector) : 

T 

H2  =  \  II  +    f*(f  '  £'X  -}  '  (3.7.37) 

The  necessary  conditions  for  optimization  are 

T 
H.    =  0* 

*'*        ~ 

and 

T 
e*  =  ~  (Ho   > 

or,  after  performing  the  differentiations 

Hn   =  el   +  er=  0T  (3.7.38) 

'  e 

and 

f*  =  ?'X  ^   •  ■  (3.7.39) 

The  boundary  conditions  are 

e  (0)  =  e  rtj  =  0  (3.7.40) 

_*       _*   f     _ 

because   \(t)  is  free  at  t =  0   and  at   t  =  t  . 

Using  Equation  (3.7.38)  in  (3.7.39)  and  (3.7.40),  one  obtains 

e  =  f ,v  e  (3.7.41) 

and 

e(0)  =  e(t  )  =  0  .  (3.7.42) 

Clearly,  Equations  (3.7.36),  (3.7.41)  and  (3.7.42)  form  a  linear  TPBVP 
in  X   and  e.   This  problem  should  be  easier  to  solve,  in  theory,  than  the 
original  TPBVP  of  the  state  and  adjoint  variables  and  should  give  an 
optimal  choice  for  the  multiples  of  X. 
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In  several  attempts  this  linear  TPBVP  could  not  be  solved  for 
the  guessed  state  trajectory  by  either  the  quasilinearization  program  or 
by  the  IBM  scientific  subroutine  package  program  DLBVP. 

The  prime  reason  for  the  failure  of  the  method  of  quasilineariza- 
tion was  finally  found  to  be  the  numerical  inaccuracies  in  the  compu- 
tations which  were  dominant  in  spite  of  using  the  double  precision 
arithmetic.   The  problem  was  perpetuated  and  amplified  by  the  large 
numbers  in  the  right  side  of  the  state  and  adjoint  equations  and  in  the 
derivatives  of  these  equations.   The  matrices  to  be  inverted  at  the 
various  stages,  one  at  every  step  of  integration  and  another  at  the  end 
of  the  integration,  were  ill-conditioned  for  inversion.   When  a  dif- 
ferent subroutine  for  inversion  than  what  came  with  the  QUASI  program  was 

m 

used,  different  numbers  resulted.   The  entries  of  the  matrices  T  and  Tf   D. 

i=l  * 

were  very  large  and  resulted  in  large  numbers  for  some  entries  of  C. 
This  made  the  matrix  C  ill  conditioned  for  inversion.   Any  error  in 
inverting  the  matrix  C  would  be  amplified  in  the  values  of  e    . 
This  amplification  was  due  to  the  multiplication  of  the  inverse  of  C 

by  T,  a  matrix  used  in  generating  e  ..  .   If  e    was  in  error,  all  other 

J  -'  b  -m+1      -m+1 

e  would  be  in  error  because  they  were  generated  from  e    . 
-i  -m+1 

If  the  inversion  of  C  was  accurate,  then  the  first  six  entires 

of  e   computed  from  Equations  (3.7.16,17,18)  should  be  almost  zero.   But, 

3 
instead,  large  values  of  the  order  of  10  were  obtained!   This  obviously 

indicated  that  ev  and  hence  all  the  e.  were  being  computed  inaccurately. 
-1  -i 

The  quasilinearization  program  failed  to  solve  the  human  motion 
problem  due  to  the  above  three  reasons,  and  primarily  due  to  the  last 
one,  the  numerical  inaccuracies.   This  was  felt  to  be  rather  difficult 
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to  overcome  since  it  was  intimately  related  to  the  method  used  to 
solve  the  linear  TPBVP  and  the  structure  of  the  original  TPBVP.   So 
far  as  the  method  was  concerned,  the  key  point  was  that  the  matrix  C  was 
becoming  ill  conditioned.   This  occurred  because  the  recursive  rela- 
tionship between  e.  and  e.  .    has  been  used  to  generate  a  relationship 
r         -l     -l+l 

m 

between  e   and  e  .  ,  which  resulted  in  the  product  Tt  D  with  large 
-1     -m+1  i=l   x 

m 

entries.   A  look  at  the  expression  for  C  in  terms  of    Tf  Di ,  B± , 

and  B   (Equation  (3.7.8))  would  make  it  clear  that  with  such  a  value 
r 
m 
of     TT  D  ,      C  would  automatically  be  ill  conditioned. 

i=l  ± 

The  other  standard  methods  of  solving  linear  TPBVP  (Equations  (3.7.5), 

(3.7.5)),    for  example,  the  transition  matrix  algorithm  might  have  been 

numerically  more  stable.   However,  other  difficulties  arose  because  these 

methods  require  several  forward  integrations  in  one  iteration.   This 

means  calling  the  subroutine  DIFEQ  (the  subroutine  to  generate  the  right 

side  of  the  differential  equations  and  its  derivatives)  many  more  times. 

This  increases  the  computing  time  enormously. 

With  a  step  size  small  enough  for  the  integration  to  be  stable, 
the  storage  requirements,  computing  time,  and,  therefore,  the  cost  of 
computing  increases  considerably.   Even  if  these  factors  are  absorbed, 
it  may  still  be  necessary  to  try  several  initial  guesses  of  the  adjoint 
variables  to  get  the  method  to  converge. 

In  view  of  the  above  problems,  it  was  concluded  that  the 
standard  method  of  quasilinearization  was  not  suitable  for  the  human 
motion  problem  and  so  should  not  be  pursued  further. 
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3.8.   Steepest  Descent  Methods  for  Solving 
the  Minimum -Time  Kip-Up  Problem 

Three  different  formulations  of  first-order  steepest  descent 
methods  were  used  after  the  method  of  quasilinearization  was  unsuccess- 
ful in  solving  the  minimum-time  kip-up  problem.   These  formulations 
differ  from  each  other  in  the  construction  of  the  cost  functional, 
handling  of  the  terminal  constraints,  treatment  of  the  control  incre- 
ments and  in  the  method  of  adjusting  the  final  time.   The  basic  features 
of  these  three  formulations  are  described  below. 

Formulation  1 

a.  The  cost  functional  is  the  final  time.   The  terminal  errors 
and  the  cost  functional  are  reduced  simultaneously. 

b.  The  adjustment  of  the  final  time  is  done  by  extending  or 
truncating  the  final  end  of  the  trajectories. 

c.  The  control  functions  take  the  form  of  a  sequence  of  constrained 
and  unconstrained  arcs.   Improvements  are  made  at  the  uncon- 
strained parts  only,  including  the  junctions  of  the  constrained 
and  unconstrained  arcs. 

The  method  is  based  on  the  works  of  Bryson  and  Denham  [27,28] 
and  Bryson  and  Ho  [32]. 

Formulation  2 

a.  The  cost  functional  is  the  sum  of  a  scalar  representing  the 
final  time  and  a  norm  of  the  terminal  error. 

b.  A  change  in  the  independent  variable  t  is  introduced  by  defin- 
ing the  transformation 
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t  =  ai 
where  a  is    a   constant,    or 

**  =  o 

dT 
so  that  a   is  treated  as  an  additional  state  variable.   The  final 
time  t   is  directly  proportional  to  a   when  T   is  held  fixed. 
Long  [29]  used  this  transformation  of  the  independent  variable 
to  convert  free  end  point  TPBVP  to  fixed  end  point  TPBVP  for 
solution  by  the  method  of  quasilinearization. 
The  cost  functional  is  reconstructed  as 

2     6  2 

J  =  K  a  +   £  K.§.  ,      K  ,K.  ,.  .  .  ,K  >  0. 

o      .  H  l  l        o   l       6 
1=1 

No  terminal  constraints  were  introduced  in  this  formulation 
since  they  (the  § .  ' s)  are  included  in  the  cost  functional. 
c.   The  control  functions  are  assumed  to  be  free  to  change  in  any 
direction  while  computing  the  gradient  of  the  cost  function. 
That  is,  the  control  constraints  were  not  considered  when  com- 
puting the  gradient  of  the  cost  functional  used  to  find  a  suit- 
able increment  in  the  control  and  a.      The  control  constraints 
were  imposed  in  the  next  iteration  during  forward  integration 
of  the  state  equations.   When  the  computed  control  violated 
a  constraint  in  a  subarc  it  was  set  equal  to  its  limit,  the 
constraint  function  on  the  subarc.   This  approach  for  treating 
constraints  on  the  control  has  been  used  by  Wong,  Dressier  and 
Luenburger  [33], 
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Formulation  3 

This  formulation  consisted  of  the  features  (a)  and  (b)  of 
Formulation  2  and  (c)  of  Formulation  1. 

The  derivations  of  the  numerical  algorithms  for  these  formula- 
tions are  now  presented.   The  basic  concepts  on  which  these  algorithms 
are  based  are  available  in  the  literature  [28,32].   The  results  are 
derived  in  a  manner  suitable  for  analysis  of  the  outcome  of  the  numer- 
ical computations,  and,  thus,  the  derivations  presented  here  are  slightly- 
different  from  those  found  in  the  literature  for  any  gradient  method 
approach  to  the  computation  of  optimally  controlled  motion. 

3.8.1.   Derivations  for  Formulation  1 

Suppose  we  have  a  continuous  nominal  control  u  (t)  and  u  (t) 
and  a  nominal  final  time  t  .   These  control  histories  have  some  parts 
lying  on  the  constraints  SJ (X) ,  i  =  1  or  2,  and  the  remaining  parts  lie 
away  from  the  constraints.   The  parts  lying  on  constraints  will  be 
called  "constrained  arcs"  and  the  parts  lying  off  the  constraints  will 
be  called  "unconstrained  arcs."   The  interesections  of  constrained  and 
unconstrained  arcs  will  be  called  "corner  points."   A  nominal  guess  for 
the  control  variables  consists  of  specifying  the  corner  points   and 
values  of  the  control  at  unconstrained  portions.   On  the  constrained 
arcs  control  variables  are  generated  from  constraint  functions.   An 
initial  choice  of  the  control  history  and  the  final  time  will  not,  in 
general,  satisfy  the  boundary  condition  and  will  not  do  so  in  minimum 
time  either.   One  can  improve  upon  the  trajectories  in  the  following  way. 

At  first,  we  establish  how  a  particular  state  variable  X   (the 
i   component  of  the  state  vector  X)  changes  at  the  final  time  due  to 
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a  small  change  in  the  control  history  and  the  final  time.   For  that 
purpose,  let  a  cost  functional  be  defined 

J  =  Xi(tf)  •  (3-8-1) 

Let  \X(t)  be  an  arbitrary  time-varying  vector  of  dimension  six. 

Since  the  system  satisfies  Equation  (3.1.1),  the  final  value  of  the 

state  variable  will  also  be  given  by 

t 

j'  =  X^tJ  +  J    X1   (f-X)  dt  .  (3.8.2) 

f    "'o   -    -  - 

If  the  control  variables  u  (t)  and  u2(t)  change  by  a  small 
amount  6U  (t)  and  6u  (t) ,  there  will  also  be  a  small  change  in  the 
state  variable  X(t) ,  denoted  by  6X(t) ,  throughout  the  trajectory. 
It  is  clear  that  these  changes  in  the  control,  denoted  by  6u(t)  and 
in  the  state  variables  denoted  by  6x(t) ,  will  not  be  independent  of 
each  other.   Apart  from  changing  the  control  history,  an  increment  to 
the  final  time  t   by  a  small  amount  St   and  small  increment  6x(0)  to  the 
initial  state  vector  are  also  prescribed.   The  first-order  change  in  J 
due  to  the  changes  in  the  control  and  the  final  time  is  given  by 

AJ'  =  X\tf)  dtf  +  {ex1  -  (\1)T  6x}t  +  (XV  fiX(O) 

+  /  f  {(^>T  !»x  +  (^}}  6-  dt  +  J*      (-±)T  !'u  6-  dt' 

(3.8.3) 
The  Su  is  chosen  in  the  following  way:   For  the  unconstrained  parts, 
6u  is  completely  free.   The  parts  presently  on  the  constraints  will 
remain  on  the  constraints  for  the  same  periods  of  time  as  before. 
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For  these  portions,  the  change  in  control  is  given  by  the  shift  of  the 
constraint  due  to  state  changes  according  to  the  relation 

6u  (t)  =  S3         4x(t)  .  (3.8.4) 

Let  C  denote  all  those  portions  of  the  trajectory  of  the 
i 

control  u.(t),  i  =  1  or  2,  which  lie  on  any  of  the  constraints  S.  or 
S..   Also,  let  C.  denote  all  those  portions  of  the  control  u.  which 
do  not  lie  on  a  constraint.   If  the  expressions  for  $u  on  the  con- 
strained arcs  given  by  Equations  (3.8.4)  are  substituted  in  Equa- 
tion (3.8.3)  and  the  integration  of  the  last  term  of  the  right  side  of 

o   1 
Equation  (3.8.2)  is  split  into  integrations  over  the  intervals  CL  ,  C^, 

C  ,  and  C  ,  one  obtains 

6j'  =  xi(t^)dt„  +  {ex1  -  aV  &x},_,    +  {a1)7  5x} 

1   x   *■     "     J  X~Xf  t=o 

tf 

+  Jo    {(^1)T  !'X+  ^l  5-  dt 

+  /  (xV  f,    6u   dt  +  /  (X1)1  f ,    S3        6x  dt 

r°  X         r1  X  - 

Ll  Ll 

+  P   (\V  f .  ,   6u  dt  +  J    (X1)1  t  ,n      4        OX  dt  .  (3.8.5) 

o   "     "   2    2      1  "    "   2    'X 
C2  C2 

If  X  (t)  is  computed  such  that 

\*(t,)  =  1      for   i  =  j  (3.8.6) 

=  0      for   i^  j 
where    X.  =  j    element  of  X 


and 


-="^*=iS\tS!S'!J'  -1        <3'8,7) 


SI 


where 


6  =  0  on  C°     §   =  1  on  C„ 
1  1'     1  1 

S2  =  0   on  Cg.     *2  =  1   on  C* 


and 


u1    and  u    used  in  Equation  (3.8.7)  are  computed  only  when 

I  I 

the  controls  u,  and  u   are  on  a  constraint  and  u   and  u   are  replaced  by 

the  constraint  expressions  S.(X),  one  obtains 

fij'  =   x\t  )  dt  +  (X1(0))T  6x(0)  +  /.(X1)1  f,„   K  dt 
11  Cl         1 


i.T 


+  Jo(^   *•«   S  dt 


=  f*(X  (t  ),u(t  ))dt  +  (^(O))1  6X(0)  +  J*  (XV  f,    «ux  dt 

cx       1 

+  .f0(^)T  *,u   *u2  dt  <3-8"8) 

2 

where  f   is  the  i    element  of  the  vector  f(X,u). 

Equation  (3.8.8)  is  the  desired  expression  for  the  change  of 
the  state  variable   X1  at  the  final  time  t   due  to  (1)  a  small  arbitrary 
increment  §u(t)  given  to  the  control  variable  u(t)  over  the  unconstrained 
portions,  (2)  a  small  increase  dt   of  the  final  time  tf ,  and  (3)  an 
arbitrary  small  change  $X(0)  of  the  initial  state  vector  X  . 

Similarly,  one  can  find  the  expressions  for  the  change  in  the 
terminal  value  of  any  other  state  variable  X  .  It  can  be  seen  that  if 
one  constructs  the  (6x6)  matrix 
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R(t)  ^[^"(t),  X2(t),  \3(t),  \4(t),  X5(t),    X6(t)J 


so  that  R(t)  satisfies 


R(t  )  =1      (6x6  unit  matrix)  (3.8.9) 


and 

R(t) 


-  [*>x  +  6l£'Ul  Ul,x+  62  f->u2   u2,x  ]   R(t)        (3-8-10> 


where  the  meaning  of  the  various  terms  in  the  parentheses  of  right 
side  of  Equation  (3.8.10)  is  the  same  as  that  in  the  same  terms  appear- 
in  Equation  (3.8.7),  then  the  change  in  the  terminal  value  of  the  state 
vector  is  given  by 

§X(t  )=f  dt_  +  RT(0)OX(0)  -f   RT  f,   Su  dt  +  f   RT  f,    §uo  dt  . 

-  f    -   f        -     *•  o    -  u_   1    u  o    -  u„    2 

Cx        1        C2       2 

(3.8.11a) 
If   we   choose  that   6X(0)  =  0,  we  obtain 

6x(tJ  =  f  dt  +  f  RT  f,    5u  dt+f  RT  f,    6u   dt.  (3.8.11b) 

-  f    -   f  Jco    -  ut        1    J,o    -  u2   2 

Following  the  method  prescribed  by  Bryson  and  Ho  [32] ,  it  will 
now  be  attempted  to  make  improvements  in  the  terminal  errors  given  by 

|  =  X(tf)  -  Xf 

and  at  the  same  time  reduce  the  final  time  t  .   Thus,  since  it  is  being 
sought  to  minimize  t  ,  one  would  maximize  -dt  or  minimize  dt  with 
respect  to  ^u   6u   subject  to  the  constraint 

6x(t  )  =  f(t  )dt   +  f   RT  i,    6u   dt+f   RT  f,    $uQ  dt       (3.8.12) 
-f    -  f   f   **  o    -  u,    1      °  o    -  u„   2 

Cl        X  C2        2 


cSiJ 


with 

°X(tf)  =  6xP(tf) 

.  p 
where  oX  (t  )  is  a  chosen  decrement  in  the  terminal  error  such  that 

&u  maintains  the  first-order  approximations. 

In  this  incremental  minimization  problem,  the  incremental  cost 

functional,  dt  ,  to  be  minimized,  and  the  constraints  are  linear  in  the 

incremental  control  parameter  6u.   Such  a  problem  does  not  have  an 

extremum.   However,  since  these  are  linearized  relations  obtained  from 

a  nonlinear  system,  the  increments  6u  ,  &u 0>  and  dt   should  be  small  for 

X      £  X 

the  first-order  approximation  to  be  valid.   To  limit  the  increments 
&u  ,  8u  ,  and  dt  ,  the  following  quadratic  penalty  term  to  the  incre- 
mental  cost  dt   is  added: 

|  b  dtj  +  \  Jq   6uJ  Wx(t)  dt  +  |  JQ  ?>u22   W2(t)  dt,         (3.8.13) 

'  Cl  C2 

where  b  is  a  positive  scalar  quantity  and  W  (t)  and  W  (t)  are  positive 

X  Ci 

scalar  quantities  specified  as  functions  of  time. 

Adding  these  quantities  to  the  cost  functional  and  adjoining 
the  constraint  relations  (3.8.12)  to  the  resulting  expression  by  a 
multiplier  v  (a  six-dimensional  vector) ,  the  following  problem  is 
obtained. 

Minimize  wrt,  &u  ,  6u  ,  and  dt 

X       £1  X 

k4bdtf4i0  s2  widt4/o du'  w2dt+vT{f(tf)dtf 

Cl  C2 


RTf,         §u1dt+   f      RTf,         6U   dt    -   cXP(t„)h 
o        -  u„         1  i,o        -  u„        2  -        f  J 


+  1 


(3.8.14) 
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If  the  derivatives  of  this  functional  with  respect  to  §u  ,  ^u9,  and 


dt   are  set  to  zero,  one  obtains 
f 

dt     =    -  i   [l+vT   f(t    )]    , 
f  b  —     -      I 


(3.8.15) 


S  -  -  ^  ■  'T^  rv 


on     C 


°       ^ 


and 


--i-fT 

U2  = 

Rv  on      Cr 


J 


(3.8.16) 


for    an  extremal. 

Using   Equations    (3.8.15)    and    (3.8.16)    in    (3.8.12),    one  obtains 


6XP  =    -  I     f{l  +  vTf} 


-   I.  i    v 


(3.8.17) 


or 


where 


v   =   - 


1  T 

I.     +-  f   f'(t„) 

ijn{f     b  -   -        f 


-l  r 


•r<V  +  b*(V 


(3.8.18) 


I.,  =  r     RTf,      W"1  fT     Rdt+f     RT  f ,      W       fT        R  dt. 
M     Jo        -  ux    1      -'ux  J,o  -'u2   2      -  u2 

(3.8.19) 

The  value  of   &XP(t    )  ,    the   desired   change   in  the   terminal   values  of   the 
state  variables,  may  be   chosen   as    a  decrement    in   the   terminal    error  §. 

6xP(tJ    =    -    e.    §.  (3.8.20) 

if  11 

where 


and 


0  <  e.   £  1 

i 


§XP(t    )    is    the    i        element   of   6xP. 
if 
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The  steepest  descent  algorithm  for  Formulation  1  can  now  be 
described  as  follows: 

1.  Guess  a  nominal  control  history  u  (t)  and  u  (t)  and  a  final 
time  t  . 

2.  Integrate  the  Stcte  equations  forward  with  the  nominal  control 
and  the  nominal  final  time  with  the  initial  values  of  the  state 
vector  given  by  X(0)  =  X  .   Store  X  (t) .   Compute  and  store  |. 
Find  the  norm 

6      ,03 


Hill  -U .  k/J 


i=l 
for  some  positive  K.,  i  =  1,...,6  (to  be  specified). 

3.  Save  the  controls  u   and  u   in  another  variable  L'OLD,  the  corner 

J-  Ci 

times  N  in  the  variable  NOLD  and  the  final  time  t   in  TFOLD. 

4.  Set  R(t  )  =  I,  the  (6  x6)  unit  matrix.   Integrate  backward 

T 
Equations  (3.8.10).   At  the  same  time  compute  and  store  R  f, 

Ul 

T 
and  R  f,   .It  is  not  necessary  to  store  R(t) . 

"  U2 

5.  Select  the  positive  constant  b  and  the  positive  quantities 

W1 (t)  and  W  (t).   Unless  there  is  some  special  reason,  W  (t) 
and  W  (t)  may  be  taken  as  a  positive  constant  equal  to  W  (to  be 
specified).   In  such  a  case,  the  storage  for  W  (t)  and  W  (t) 
will  not  be  needed.   Select  e.,  i  =  1,2,. ..,6  such  that 

0  <  e  £  1. 

1 

6.  Compute  I , ,  using  Equation  (3.8.19). 

7.  Compute  v,  6u  ,  §u    ,    and  dt   from  Equations  (3.8.18),  (3.8.15), 
and  (3.8.16).   Note  that  Su   and  su   are  defined  for  the 
unconstrained  parts  C   and  C   only,  and  are  to  be  computed  only 
on  those  parts. 
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Change  the  final  time  tf  to  tf  +  dtf.   If  the  integration  step 
size  is  h,  this  means  that  the  total  interval  is  to  be  increased  or 
decreased  at  the  final  end  (depending  upon  whether  dtf  is  positive  or 
negative)  by  dtf/h,  rounded  to  the  nearest  integer  number  of  integra- 
tion steps.   m  the  case  of  a  positive  dtf,  the  controls  u^t)  and 
u2(t)  in  the  interval  TFOLD  *  t  *  tf  (new)  will  be  given  by  extrapola- 
tion of  the  existing  curves.   if  a  control  function  is  on  a  constraint 
at  the  final  time,  TFOLD,  it  should  remain  on  the  same  constraint  in 
this  period  and  is  to  be  generated  during  the  forward  integration  in  the 
next  steps. 

8.   Set  the  control   u.  =  UOLD.  +  6U.,  i  =  i  or  2,  for  the  uncon_ 
strained  portions.   Integrate  the  state  equations  forward  with 
the  new  control.   The  corner  times  N  are  to  be  generated  again 
during  the  forward  integration  as  described  graphically  in 
Figure  11. 

Figure  11  shows  the  different  cases  which  may  arise  during  the 
forward  integration.   It  is  seen  that  a  new  corner  time  is  generated  at 
the  point  where  the  unconstrained  u.  +  5U.  curve,  extrapolated  if  neces- 
sary, meets  the  constraint.   This  part  of  the  forward  integration. makes 
the  programming  complicated  with  many  logical  program  statements. 

9.   Find  the  errors  in  the  boundary  conditions  §  and  find  the  norm 
I!  Sjl     as  in  step  2.   If  this  norm  is  less  than  that  of  the 
previous  iteration  go  to  step  3  to  continue  the  iterations. 
If  it  is  not  reduced,  then  do  the  following: 
a.   If  dtf  is  too  large,  increase  b,  and  go  to  step  9c  bel< 
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What  is  "too  large"  is  to  be  found  by  trials.   In  the 
present  case  more  than  five  integrations  steps  were 
found  to  be  too  large  for  a  value  of  dt  . 

b.  If  dt   is  not  too  large,  decrease  the  e.  and  go  to 
step  9c  below. 

c.  Set  u  =  UOLD,  t  =  TFOLD,  N  =  NOLD  and  go  to  step  5   above. 
The  stopping  conditions  are  given  by: 

1.  ||  §J|  is  a  small  number,  meaning  that  the  boundary  conditions 

have  been  adequately  met. 

T 

2.  (1 + v  f)    is  negative  or  zero.   This  means  that  the  correc- 

tf 
tions  to  the  final  time  which  would  be  computed  in  the  next 

iteration,  using  Equation  (3.8.15)  would  increase  the  final  time. 

T  T 

3.  For  i  =  1  and  2  the  quantities  v  R  f,   are 

l"  Ui 

a.  positive  when  u.  is  given  by  S.  meaning  that  only  a  viola- 
tion of  the  constraint  will  improve  the  cost 

2 

b.  negative  when  u.  is  given  by  S. 

c.  Zero,  or  a  small  number  when  u.  is  not  on  any  constraint,  so 
that  6u.  computed  by  Equation  (3.8.16)  should  be  zero  or 

an  acceptably  small  number. 

Comments  on  the  Algorithm 

The  forward  integration  of  the  state  equations  at  every  iteration 
requires  special  treatment  because  of  the  presence  of  the  constrained 
and  unconstrained  arcs.   The  integration  formulas  are  slightly  differ- 
ent for  cases  when  the  right  sides  of  the  differential  equations  are 
time  varying.   On  an  unconstrained  portion  u  (t)  and  u  (t)  are 
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prescribed  as  time  varying  functions,  but  on  a  constraint,  they  are 
expressed  as  functions  of  state  variables.   This  difference  is  to  be 
taken  into  account  in  the  integration  scheme. 

The  control  constraining  boundaries  are  not  known  before  the 
state  trajectory  is  known.   Since  the  state  trajectory  is  determined 
by  forward  integration,  the  determination  of  the  corner  points  is  dif- 
ficult after  the  control  is  changed  at  of f-the-constraint  portions. 
For  this  purpose  a  first-order  prediction  of  the  location  of  the  con- 
straint and  an  extrapolation  of  the  unconstrained  u  +  6u  curve  was  made 
to  determine  if  they  intersect  to  form  a  corner  point  one  step  ahead  in 
the  forward,  state  equation  integration.   This  was  needed  because  the 
integration  scheme  used  four  consecutive  stored  states  to  generate  the 
next  state.   Much  more  computing  time  would  be  required  by  going  back 
to  an  earlier  point  and  redoing  the  integrations.   First-order  predic- 
tions of  the  states  X„  and  X_  for  the  step  K+l  were  made  from  the  step 

o      5 

K  by  the  formulas 


and 


K+l    YK      K 

x3   =  x3  +  h  x*4 


X5+1  "  A  +    h   X6  • 


h  =  integration  step  size  (constant). 

Only  X_  and  X_  are  needed  to  find  the  constraints  S. . 
3      5  l 

It  is  to  be  noticed  that  the  modification  of  a  corner  point  as 
illustrated  in  Figure  11  has  not  been  accounted  for  in  the  derivation 
of  the  algorithm.   The  errors  in  the  expression  of  6x(t  )  caused  by 
this  derivation  in  the  algorithm  from  the  theoretical  derivation  will 
be  small  only  if  the  constrained  and  unconstrained  arcs  are  much  longer 
than  the  shift  of  the  corner  point.   This  is  a  significant  source  of 
error  when  the  unconstrained  arcs  are  small. 
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3.8.2.   Derivations  for  Formulation  2 

In  this  formulation,  we  transform  the  independent  variable 

t  =  0T  (3.8.21) 

where  T  is  the  new  independent  variable. 

Then  the  system  equations  become 
dX 


_  =  x'  =  a   f  (X,u) 

dT   -      -  -  - 


da    '   „ 

—-=   a     =0 

dT 


(3.8.22) 


where  (')  denotes  differentiation  with  respect  to  T.  (3.8.23) 

Boundary  conditions  are 

X(0)  =  X     (given) 


with 


X(T  )  not  specified  and  to  be  determined, 


and 

t  ,  the  fixed  final  value  of  T,  given. 

The  cost  functional  to  be  minimized  is 

J  =  -  £  K  §2  +  K  a2  (3.8.25) 

2  .  _   11     o 
1=1 

with  the  values  of  the  K   and  K.  specified.   Constraints  on  the  control 

o       1 

given  by  Equation  (3.1.7)  remain  unchanged. 

The  steepest  descent  algorithm  for  minimizing  the  cost  functional 
of  Equation  (3.8.5)  was  constructed  as  follows. 

Let  X(t),  a  vector  of  dimension  6,  and  T](T),  a  scalar,  be  arbi- 
trary functions  of  T.   Since  the  system  satisfies  Equations  (3.8.22), 
the  cost  functional  is  also  equal  to 


9! 


6    i-  -i  2  f  t 

j'=i  ZK.   Xi(tf)-X^J  +KQa2  +  /   >T[cvf  -x']dT-  /   V  dr.   (3.8.26) 


In  this  expression  §   has  been  replaced  by  (X  (tf)  -  Xf )  . 
Integration  by  parts  yields 


J  =  9  + 


where 


/    '   X.V  dr  -  [XTX]Qf  +  jQf  \'X  dt  -  [V]0f  +  /Qf  71'   adx 


1  i  i   2        1  2 

»-£  E  Ki[x  (V-V  +  2  V  • 

i=:l 

For  first-order  variations  in  u  and  X 

t. 


+  cp,      §a  +    f        \T[f   §o/  +   a  f,„    6X  +   o  f .      flu]    dT 


-  CXT  6X]       +  J       X/T  Ox  di 
"  rf         o 


[7]       _T)    ]    Scv  +    f  tj'    Oa  dT.       (3.8.27) 
t  „        o  u 


Let    \(T)    and  T)(T)    be  chosen  such   that 


X    (Tf)    =  cp,x 


X'(T)   =   -   (a  f,x)      X 


71' (T)   =   -   XXf 


H(V  =cp,a 


(3.8.28a) 

(3.8.28b) 

(3.8.28c) 
(3.8.28d) 


Then, 


6j 


'  -  T)     fla  +  J       XTa  f ,      $u 


dT 


(3.8.29a) 


6a  +   a  f        X      [f ,        5u     +   f ,        flu   ]    dT, 
o  J         -        -  u„         1        -  u0        2 


(3.8.29b) 
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Equations  (3.8.29)  give  the  gradient  of  the  cost  function 
(3.8.26)  with  respect  to  a   and  6u(T).   Let  us  choose 

ba  =  -  e  T]   s  (3.8.30) 

a  o 

6U;L(t)  =  -s  \T  f,u   Sl(T)  (3.8.31a) 

6u(t)  =  -s  XTf,     e  (t)  (3.8.3lb) 

z  -  -   u2  t 

where  e  ,  e  (T)  ,  S-CO,  and  s  are  positive  quantities.   It  can  be  seen 
a   1      2 

that  if  these  quantities  are  chosen  such  that  6cv,  6u  (t)  ,  and  6u  (T)  are 
sufficiently  small,  then  the  first-order  approximations  are  valid  and  the 
reduction  of  the  cost  is  guaranteed.   The  purpose  of  making  e^(T)  and 
e  (t)  functions  of  T  is  to  have  control  over  the  amount  of  improvements 
§u   and  6u   during  the  time  interval.   They  will  normally  be  1.0  all 
over,  except  for  special  reasons.   The  values  of  e  ,  e  (T) ,  and  e2(T) 
are  to  be  specified  independently.   Once  they  are  decided  upon,  s  may 
be  calculated  on  the  basis  of  how  much  improvement  to  the  cost  is 
desired  in  one  iteration  while  remaining  within  the  linear  approxima- 
tion's valid  range. 

Suppose,  in  any  iteration,  the  cost  is  J,  and  it  is  desired  to 
reduce  the  cost  by  e J  in  the  next  iteration,  where  0  <  s  S  1.   With 
the  choice  of  6u  (T)  ,  5u  (T)  ,  and  §cv  already  made  above,  Equation 
(3.8.29b)  yields 

-ej=-s  e  if-saf   [e  (T)(^Tf,   )2  +  e  (T)(\Tf,   )2]  dT.    (3.8.32) 
ao     J  1     -  -  U„       2     --u^ 

o  1  2 
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Defining 

f 


=  a  f   [e1(T)(XTf,   )   +  e9(T)(\Tf   )-]  dT  (3.8.33) 

«J        1      —  —  U,         A  —    —    U 


B 

1  2 


one  obtains 


-ej  =  -s[e  t\2   +  B]  (3.8.34) 

a  o 

or 

s  =  ej/(s   if  +  B).  (3.8.35) 

a  o 

Equation  (3.8.35)  gives  the  value  of  s  required  in  the  determination 

of  6^  (T)  ,  6u  (T) ,  and  6a  via  the  gradient  terms,  which,  under  first- 

order  predictions,  would  reduce  the  cost  by  ej.   In  actual  computation, 

a  different  change  in  the  cost  is  to  be  expected  because  e  will  be 

finite  and  the  first-order  prediction  will  not  be  exact.   This  is  true 

because  the  system  is  nonlinear,  and,  because  in  some  intervals  of  T, 

6u  may  not  be  implemented  due  to  the  constraints. 

If  the  cost  functional  is  not  improved,  then  e  is  to  be  reduced 

and  6u  and  6q<  are  to  be  recomputed.   If  the  cost  function  is  reduced 

but  the  terminal  error  is  increased  (meaning  thereby  that  the  cost  has 

decreased  by  reducing  a   but  increasing  the  terminal  error) ,  K   is  to  be 

reduced. 

The  selection  of  e  was  done  in  the  following  way.   The  change 

a 

in  the  cost  due  to  the  change  in  a   and  the  change  in  u(T)  are  given  by 

-s  s  T)   and  -sB,  respectively, 
a  o 

Let 

-s  e  tl   =  e  (-sB)    where  e   :>  0  (3.8.36) 

a   o     a  a 

so  that  the  change  in  cost  due  to  change  in  a   is  e   times  that  due  to 

a. 

the  change  in  u(t).   Then  e   is  given  by 
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e  =  e   B/T|  (3.8.37) 

a    a    o 

and  then  from  Equation  (3.8.34),  s  will  be  given  by 

s  =  ej/[B(l+  e  T,  )}  .  (3.8.38) 

a  o 

Since  it  is  intuitively  easier  to  specify  €   than  e  ,  the 

advantage  of  choosing  e   in  this  way  is  obvious. 

3.  , 

The  Stopping  Conditions 

The  exact  stopping  conditions  for  the  computations  should  also 
correspond  to  the  conditions  for  optimality  and  they  are  suggested  by 
Equations  (3.8.29).   Clearly,  when  the  following  three  conditions  are 
satisfied: 

1.  T)  =0   (or  small)  , 

o 

T  T 

2.  X  f,    and  X  f ,    are  zero  (or  small)  when  the  control  variables 

"  -  Ul      "  "  U2 

u   and  u  ,  respectively,  are  not  on  a  constraint, 

T         T 

3.  when  u  or  u   are  on  a  constraint,  X  f ,    or  X  f,   ,  resoec- 

1     2  -  -  u     -  -  u2 

tively,  are  of  opposite  signs  to  those  of  the  constraint,  i.e., 
positive  on  a  lower  constraint  and  negative  on  an  upper  constraint; 
then  no  improvement  on  cost  is  possible  without  violating  the  constraints. 
These  indicate  stopping  condition.   While  performing  the  last  check  for 
stopping,  the  adjoint  equations  must  be  modified  according  to  Equation 
(3.8.41)  of  Formulation  3  when  any  portion  of  the  control  trajectories 
lies  on  a  constraint.   This  is  necessary  in  order  to  take  into  account 
the  dependence  of  the  control  constraints  on  the  state  variables. 
However,  there  are  several  problems  to  expect  from  these 
stopping  conditions.   These  conditions  will  not  be  satisfied  even  if 
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a  trajectory  is  slightly  different  from  the  optimal.   The  adjoint 
variables  are  extremely  sensitive  with  respect  to  changes  in  the  control 
functions  and  a  slight  imperfection  in  the  solution  will  yield  a  signif- 
icantly different  set  of  adjoint  variables.   They  will  not  meet  the 
stopping  conditions  described  above.   Also,  the  first-order  methods  show 
poor  convergence  near  the  optimal  solution  and  may  never  go  sufficiently 
close  to  the  exact  solution  in  a  reasonable  amount  of  computation. 
Added  to  these  problems  is  the  problem  of  computation  errors  (round  off 
and  truncation  errors)  which  grow  in  integrations  over  a  lengthy  inter- 
val, and  will  begin  to  show  up  then  the  terminal  errors  and  the  gradient 
of  the  cost  functional  are  reduced  to  small  values.   In  view  of  these 
factors  it  is  very  unlikely  that  the  stopping  conditions  will  be  pre- 
cisely obtained.   One  should  therefore  find  some  other  stop- 
ping criteria  to  decide  when  to  stop  the  computations  and  accept  an 
approximate  solution. 

Since  it  is  desired  to  minimize  the  terminal  errors,  one  can 
test  directly  whether  they  are  small  enough  or  not.   For  time  optimality, 
it  has  been  seen  that  the  optimal  solution  is  bang-bang  with  the  excep- 
tion of  singular  arcs.   It  is  very  difficult  to  start  with  a  continuous 
nominal  control  trajectory  and  improve  to  an  exact  bang -bang  control  in 
any  part  of  the  trajectory.   So,  a  steep  control  function  arc  joining 
an  upper  and  lower  constraint  should  be  an  acceptable  approximation  of 
a  bang-bang  solution.   If  there  are  singular  arcs,  they  must  be  tested 
for  optimality  by  seeing  if  they  are  stable  or  not.   Also,  T)   should  be 
small  or  negative.   A  negative  T]   indicates  that  the  cost  functional 
can  be  reduced  further  by  increasing  a.      This  corresponds  to  increasing 
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the  final  time,  which  we  wish  to  minimize.   If  we  already  have  acceptable 
terminal  errors,  no  further  increment  in  a'  is  needed,  in  this  case. 
Lastly,  with  all  the  above  conditions  satisfied,  the  final 
stopping  condition  is  the  insignificant  improvement  in  the  cost  func- 
tional by  further  computations. 

Choice  of  the  K. 

1 

If  any  K.  is  large  compared  to  the  others,  it  may  be  expected 

that  the  corresponding  boundary  error  will  be  reduced  at  the  expense  of 

the  others.   On  the  other  hand,  not  all  errors  may  improve  acceptably  if 

all  the  K   are  made  equal.   If  K   is  too  small,  the  minimization  of  a 
i  o 


w 


ill  be  slow.   If  K   is  large,  a   may  be  reduced  too  much,  resulting  in 


o 


large  boundary  condition  errors.   A  reasonable  way  to  select  Kq  is  to 


increase  or  decrease  K  with  the  current  norm  of  the  terminal  error 

o 

P1  =  Hill  =  |  £  K.  §J  .  (3.8.39) 

i=l 

The  values  of  the  K.  for  i  =  1,2,. ..,6,  may  need  to  be  adjusted  during 

the  iterations  to  achieve  faster  improvement  of  the  individual  terminal 

errors. 

The  algorithm  developed  above  can  be  summarized  as  follows: 

1.  Decide  upon  the  values  of  the  quantities  K-,  YL^,    ...  ,Kg,  e, 

e  (T) ,  e  (t)  and  e  as  discussed  above. 
1    '   2         a 

2.  Make  a  nominal  guess  of  T  ,  u^T),  ug(T)  and  a.      This  will 
require  a  few  simulation  trials. 
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3.  Integrate  the  state  equations  (3.8.22)  forward  with  the 
nominal  control  and  the  given  initial  values.   At  each  inte- 
gration step,  check  if  the  control  u  (t)  and  u  (T)  are  within 
the  allowable  bounds  (the  constraints)  which  are  to  be  computed 
at  every  integration  step.   Bring  the  controls  on  the  bounds  if 
they  exceed  the  bounds.   Store  u(T)  and  X(t) . 

4.  Calculate  the  terminal  errors  g  ,  g  , • • • ,!g  and  their  norm  P1 

and  determine  the  quantity  K  . 

1  2 

5.  Calculate   the  cost   J  =   P     +   -  K      a   . 

6.  Set   the  terminal   values   of   the   adjoint  variables: 

X    (t    )    =  K      §    ,         iHtJ    =  K  a      (from  Equations    (3. 8. 28a,  d). 
if  xi  i  o 

7.  Integrate  the  adjoint    Equations    (3. 8. 28b, c)    backwards.       Compute 

T  T 

and  store,    at   each   step,    the  values  of    X  t      CO    and   \  f         (T) . 

1  2 

8.  Calculate  B,  e  ,  and  s,  using  Equations  (3.8.33,37,38). 

9.  Calculate  6u  (T)  and  5u  (T) ,  and  &a   from  Equations  (3.8.30,31). 

JL  ^ 

10.  Update  u  (T)  ,  u  (T)  ,  and  §a   by  setting  u^T)  =  ii^T)  +  fiu^T)  , 
u  (T)  =  u  (t)  +  <3u  (T)  ,  and  a  =  a  +   &a.      Store  the  old  u  (T)  , 
u  (T) ,  and  a   in  UOLD  and  AOLD. 

11.  Integrate  the  state  equations  forward  as  in  step  3,  and  find  the 
value  of  the  cost  functional  J. 

12.  If  the  cost  functional  is  not  reduced  from  the  previous  value, 
reduce  e,  reset  the  values  of  u  (T) ,  u  (T) ,  and  a   to  the  previous 
values  UOLD  and  AOLD  that  were  stored,  go  to  step  8. 

13.  If  the  cost  is  reduced  but  P  ,  the  norm  of  the  error  has 

increased,  reduce  the  value  of  K  ,  reset  the  values  of  u  (T) , 

'  o  A 

u  (t) ,  and  a   to  the  previous  values  that  were  stored  and  go 
to  step  8. 
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14.  If  both  the  cost  function  and  P  are  reduced,  check  to  see  if 
the  stopping  conditions  have  been  met.  If  not,  go  to  step  6, 
which  is  the  next  iteration. 

15.  If  the  number  of  iterations  is  greater  than  a  specified  number, 

punch  on  cards  the  values  of  the  current  u  (T)  and  u  (T)  and 

stop.   Manually  check  to  see  if  the  program  is  performing 

well.   If  necessary,  change  the  values  of  K.,  i=l,2,...,6, 

e  ,s  ,  e  (T) ,  eo(T).   Restart  the  program  if  necessary. 
1   a   1      2 

A  Simple  Numerical  Example  of  the  Method 

A  simple  second-order  system  with  state-dependent  constraints 
on  the  control  was  taken  for  solution  by  this  method  to  examine  the 
method's  features  as  well  as  its  effectiveness.   The  system  was  defined 
as  follows: 

h   =  X2 

K=  u 

with  the  constraints  on  u  given  by 

-(0.5  +  2X2  -  0.  5X^/2  £  u  ^  (0.  5  +  2X2  -  0.  5X2>/2  . 
The  boundary  conditions  were  chosen  as 
X^O)  =  0        X2(0)  =  0 

W  = 1         X2(V  =  °- 

The  system  is  to  meet  the  above  boundary  conditions  in  minimum  time. 
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Solution 

With  the  change  of  the  independent  variable 
t  =  or 
the  state  equations  become 

Xl  =  ah 

X2  =  M 

a     =   0 
where  (')  denotes  differentiation  with  respect  to  T.   While  computing 
the  adjoint  variables  and  the  improvements  in  the  control  6u,  the  con- 
trol is  assumed  to  be  completely  free  as  mentioned  earlier  in  this 
chapter.   The  constraints  were  imposed  while  doing  the  forward  inte- 
gration with  the  improved  control  u  +  Su  in  the  next  iteration. 

The  equations  of  the  adjoint  variables  are  given  by 

T 
X  =  -X  erf  ,x 

or 

The  boundary  conditions  of  the  state  and  the  adjoint  variables  are 

^(0)  =  0     X2(0)  =  0 

X1(tf)  =  K1(X1(tf)  -  1) 

X2(tf)  =  K2(X2(tf)) 

Tl(t_)  =  K  a 
i     o 

7](0)  =  0. 
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The  cost  functional  is 

J  =  \   WV  -  I)'  +  \   K2(X2(tf))2  +  \  K/. 

In  the  attempts  to  solve  this  problem,  several  experiments  were  done 
with  respect  to  the  choices  of  the  initial  guess  of  the  final  time,  the 
quantity  K   and  e  (T). 

It  was  found  that  convergence  to  the  boundary  conditions  on 
the  state  variables  X   and  X  was  quite  rapid  from  all  the  initial 
guesses  tried,  which  may  not  be  surprising  for  this  simple  problem. 
However,  to  accelerate  minimization  of  the  final  time,  some  initial 
guesses  were  better  than  others.   It  was  found  that  if  the  final  time, 
i.e.,  the  value  of  a,    was  guessed  larger  than  optimal,  then,  in  the 
initial  iterations,  as  the  boundary  errors  are  improved  for  a  larger 
nonminimal  value  of  a,    the  control  trajectory  takes  a  shape  other 
than  optimal.   Once  the  terminal  errors  have  been  reduced  and  the  final 
time  is  still  not  minimal,  reduction  in  the  final  time  becomes  very  slow. 
Any  large  change  in  a   has  a  tendency  to  greatly  increase  the  terminal 
errors  from  their  previously  small  values. 

It  was  found  that  for  this  problem,  if  the  guessed  value  of  the 
final  time  (i.e.,  a)   was  chosen  less  than  the  actual  minimum  time,  the 
convergence  was  more  rapid.   In  such  a  case  it  was  found  that  the  cor- 
rections to  the  control  in  any  iteration  are  computed  for  a  final  time 
less  than  the  minimum  time  which  forces  the  control  to  move  to  the 
appropriate  constraint  boundaries.   At  the  same  time,  the  correction 
on  the  a   was  computed  positive,  so  that  the  correction  value  of  a   was 
approached  automatically  in  the  attempts  to  reduce  the  terminal  errors. 
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2 
Large  increases  in  a  were  checked  by  the  penalty  term  K  CK  so  that  the 

cost  functional  was  minimized  without  excessive  oscillations  in  the 

values  of  a   or  u(T). 

It  was  found  that  in  the  computation  of  the  correction  to  the 

control  by  the  formula 

6u  =  -s  e1(T)\Tfu  =  -s  ei(T)  Xg 

if  e  (T)  is  chosen  such  that  it  is  small  where  the  control  is  near 
a  constraint  boundary  and  large  when  off  the  constraint,  the  conver- 
gence to  an  optimal  bang-bang  solution  is  accelerated  after  the  control 
has  taken  a  shape  similar  to  the  optimal  control. 

Figure  12  shows  the  control  trajectories  at  different  iterations 
and  at  the  solution.   The  exact  minimum  time  solution  is  bang-bang  with 
only  one  switching  and  the  final  time  can  be  determined  by  forward  inte- 
gration of  the  state  equations  with  this  control.   In  this  way,  the 
minimum  time  was  found  to  be  3.06  seconds.   In  the  solution  obtained  by 
the  penalty  function  method,  the  terminal  errors  for  X  and  X  were  -.02 
and  -.034,  respectively,  and  the  final  time  (=a)  was  2.95.   The  values 

of  K   and  K  were  selected  to  be  1.0  for  both.   K  was  chosen  by  the 
12  o 

formula 

2     2    2 

K  =  .01  +  (g  +  §_)/ck 
o  1     z 

where  §1  =  \<Tf)  -  1>    and  ?2  =  X2(V' 

It  is  seen  that  the  final  time  in  the  solution  obtained  by  the 
penalty  function  method  is  a  little  less  than  the  actual  solution.   Such 
a  situation  is  to  be  expected  from  this  kind  of  formulation  where  a 
trade  off  between  the  terminal  error  and  the  final  time  is  bound  to 
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occur  in  order  to  reduce  the  composite  cost  function  defined  for 
this  problem. 

From  a  very  bad  initial  guess,  the  solution  was  obtained 
in  24  iterations,  taking  a  total  computing  time  of  12.88  seconds. 

3.8.3.   Derivations  for  Formulation  3 

In  this  formulation,  the  transformation  equation  of  the  inde- 
pendent variable,  the  state  equations,  cost  functional,  and  the  control 
constraint  functions  in  this  formulation  are  exactly  the  same  as  for 
Formulation  2.   Therefore,  Equations  (3.8.21)  to  (3.8.27)  of  Formulation  2 
are  also  valid  for  this  formulation.   The  difference  between  the  two 
appears  in  the  treatment  of  the  control  increments.   In  Formulation  3, 
the  control  function  is  considered  to  lie  partly  on  the  constraints,  and 
the  control  increments  are  prescribed  only  on  the  portions  of  the  con- 
trol that  do  not  lie  on  the  constraints,  as  in  Formulation  1.   The  symbols 
C  ,  C  ,  C  ,  C  for  portions  of  the  u   and  u  history  that  are  off  the 

X.  2i  X     a  X        <o 

constraints  or  on  the  constraint  will  be  used  as  in  Formulation  1. 
Considering  the  relationship  between  the  change  in  the  control  and  the 
change  in  the  state  variable  (see  Equation  (3.8.4))  as  in  Formulation  1, 
one  obtains  the  following  first-order  relationship  between  the  increment 
in  the  cost  functional  and  the  increments  in  u  and  a: 

6j'  =  T)(0)  6a  +  cv  f  XTf,   $udT+Q<  f   XTf ,    6U_  dr  .  (3.8.40) 

Jo  -  -  u   1     "  o  -  -  u    2 
Cx      1         C2       2 

The  equations  for  \(T)    and  T](t)  are  the  same  as  in  Formulation  2 
(Equation  (3.8.28)),  except  that  the  X     equation  is  now  given  by 

\'    =    ~    [f,Y+  6-,  f.    'u    +  S9  f    -u    ]T  X  .        (3.8.41) 
-         1     'X         2     'X 
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The  term  in  parentheses  on  the  right  side  of  this  equation  is 
explained  in  Formulation  1,  Equation  (3.8.7). 

Proceeding  as  in  the  derivations  in  Formulation  2,  with  obvious 
changes,  one  obtains 

B  =  «  J,,  h\u     £u   X  h   dr  +  a  /o  XT  f    fT   X«a  &       .         (3-8.42) 
C±  1     1  C2        2     2 

s  =  ^—  (3.8.43) 

C.T)  +B 
a  o 

6u,  =  -  f T   X  e  s  (3.8.44a) 

1  -  u  -  1 

6u.  =  -  f T    X  e  s  (3.8.44b) 

2  -  u2  -  2 

da  =  -  T)   s  s     .  (3.8.45) 

o  a 

The  quantities  J,  e,  e.  ,  s  ,  and  e   are  the  same  quantities  as  defined 

12       a 

in  Formulation  2. 

The  numerical  algorithm  and  the  discussions  of  stopping  condi- 
tions for  Formulation  2  and  the  comments  on  integration  with  constrained 
and  unconstrained  arcs  made  in  Formulation  1  are  valid  for  this  formulation. 

The  transformation  of  the  independent  variable  t =  qt  used  in 
Formulations  2  and  3  was  introduced  for  the  following  reason.   In  the 
usual  gradient  schemes  for  optimization  of  trajectories  and  as  in 
Formulation  1,  the  adjustment  of  the  final  time  is  done  by  extending  or 
truncating  the  final  end  of  the  trajectory.   This  was  not  considered  effi- 
cient because  the  final  time,  though  being  as  independent  as  the  control 
u  in  determining  the  optimal  trajectory,  does  not  appear  in  the  system 
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equations.  By  the  above  transformation,  a  measure  of  the  final  time,  a, 
appears  in  the  system  equations  and  therefore  direct  improvement  in  con- 
vergence was  expected. 

3.8.4.   The  Integration  Scheme  for 

the  Steepest  Descent  Methods 

The  integration  of  the  ordinary  differential  equations  for  the 
states  was  found  to  be  unstable  even  for  moderate  to  very  small  integra- 
tion step  sizes  when  the  simple  Euler  integration  scheme  was  used. 
Subsequently,  several  other  integration  schemes  such  as  a  fourth-order 
Runge-Kutta,  an  Euler  predictor  corrector,  and  an  Adams-Bashforth  pre- 
dictor corrector  scheme  were  tried.   Out  of  these  trials  the  Adams- 
Bashforth  predictor  corrector  scheme  was  found  to  be  the  best  for  our 
problem.   It  required  less  computation  per  integration  step  and  allowed 
larger  step  size  than  the  other  methods.   The  integration  step  size  was 
held  fixed  because  this  would  significantly  reduce  the  interpolation 
errors  for  the  determination  of  the  state  and  control  variables  from 
their  stored  values.   These  were  needed  for  the  backward  integration  of 
the  adjoint  equations.   The  adjoint  equations  (equations  for  the  adjoint 
variables  R,  \,    and  7])  required  a  step  size  for  integration  half  that 
required  by  the  state  equations.   The  integration  scheme  used  the  fourth- 
order  Runge-Kutta  method  for  generating  the  first  four  points  for  start- 
ing and  Adams-Bashforth  predictor-corrector  formulas  for  generating  the 
subsequent  points. 

Various  integration  step  sizes  were  tested  to  find  out  the 
largest  one  for  which  the  integrations  were  stable.   From  the  tests  the 
step  size  0.0036  second  was  found  to  be  adequate  for  the  independent 
variables  t  in  Formulation  1  and  t  in  Formulations  2  and  3  for  the 
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integration  of  the  state  equations.   It  was  also  found  that  the  integra- 
tions were  stable  if  at  the  discontinuities  of  the  right  side  of  the 
differential  equations,  the  integration  was  stopped  and  restarted  from 
that  point.   This  was  so  because  the  Adaras-Bashforth  formulas  were  not 
valid  at  a  discontinuity,  and  therefore  better  results  were  obtained  when 
starting  formulas  (four-step  Runge-Kutta)  were  used  after  the  discon- 
tinuity.  This  feature  was  therefore  introduced  into  the  program. 

3.8.5.   Initial  Guess  of  the  Control  Function 

It  was  found  necessary  for  the  initial  guess  of  the  control 
function  to  be  as  nearly  optimal  as  possible  if  considerable  computa- 
tion time  was  to  be  saved.   Simulation  (several  trial  forward  inte- 
grations) was  necessary  to  generate  an  initial  guess  for  the  relatively 
complicated  control  function.   From  the  principles  of  mechanics  it  was 
noted  that  for  a  rotating  system,  the  angular  velocity  tends  to  increase 
when  the  radius  of  gyration  of  the  system  is  reduced,  and  it  tends  to 
decrease  when  the  radius  of  gyration  is  increased.   The  gymnast  applied 
this  principle  and  counteracted  the  effect  of  gravity  by  timely  adjust- 
ments of  the  radius  of  gyration  of  the  body  about  the  horizontal  bar. 
Although  the  gymnast  is  limited  in  strength,  he  was  able  to  adequately 
shorten  or  extend  his  radius  of  gyration. 

Even  with  somewhat  of  a  physical  understanding  of  the  mechan- 
ics of  the  motion  and  the  bang-bang  principle,  it  was  not  possible  to 
devise  a  control  for  the  model  to  perform  the  desired  maneuver.   This  was 
primarily  due  to  the  fact  that  a  candidate  control  function  for  an  extre- 
mum  control  could  not  be  specified  fully  beforehand.   The  bounds  on  the 
control  were  dependent  on  the  state,  and,  thus,  were  not  known  a  priori. 
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Even  if  the  bounds  on  the  control  could  be  known  exactly,  it  is  difficult 
to  predict  the  response  of  such  a  coupled  nonlinear  system. 

The  film  data  of  the  actual  motion  were  carefully  examined  via 
the  principles  of  mechanics  but  were  of  no  further  help  in  generating  an 
initial  guess.  It  was  found  that  the  upper  limit  of  the  shoulder  torque 
for  smaller  values  of  8  was  insufficient  to  allow  the  model  to  reproduce 
the  experimentally  determined  motion.  From  the  optimality  conditions  it 
was  known  that  the  optimal  control  should  be  either  purely  bang -bang,  or 
bang-bang  with  singular  subarcs.  The  control  was  therefore  guessed  such 
that  it  coincided  with  the  constraints  on  a  major  portion  of  the  time 
interval. 

3.9.   Results  of  the  Numerical  Computations 
and  Comments 

Several  numerical  difficulties  were  faced  during  computation  by 
the  steepest  descent  methods.   These  were  caused  by: 

1.  The  system  equations  were  quite  nonlinear. 

2.  The  kip-up  motion  was  unstable. 

3.  The  control  function  had  many  switches  and  large  oscillations  oc- 
curred in  the  state  trajectory. 

4.  The  system  was  not  locally  controllable  about  several  trajectories 
that  were  tried. 

5.  Numerical  gradient  methods  have  several  problems  with  the  minimum 
time  criterion. 

The  steepest  descent  algorithms  make  control  changes  according 
to  the  local  gradient,  and  for  the  present  system  these  were  valid  only 
in  a  very  small  neighborhood  of  the  control  function.   Due  to  instability 


Ill 


of  the  perturbation  equations  about  the  kip-up  trajectory  (expressed 
quantitatively  by  large  entries  of  the  order  10  to  200,  in  the  matrix 
R(t)  in  the  earlier  segment  of  the  trajectory)  any  small  perturbation 
in  an  earlier  segment  of  the  trajectory  caused  larger  deviations  in  the 
state  variables  in  segments  near  the  terminal  time.   Thus,  unless  the 
changes  were  very  small ,  the  steepest  descent  formulas  were  not  able  to 
bring  about  the  improvement  in  the  control.   The  system  equations 
required  very  small  integration  step  size.   This  virtually  increases 
the  length  of  the  total  time  interval.   It  required  about  570  integra- 
tion steps  for  the  forward  integration  of  the  state  equations  and  twice 
as  many  steps  for  the  reverse  integraion  of  the  adjoint  equations.   With 
such  a  long  integration  interval,  the  numerical  values  of  the  gradients 
are  expected  to  be  less  accurate  due  to  accumulation  of  truncation  and 
round-off  errors.  Such  inaccuracies  near  the  earlier  parts  of  the  tra- 
jectory are  especially  detrimental  when  small  changes  in  the  terminal 
errors  are  desired  and  when  the  gradients  themselves  are  small. 

Compared  to  the  present  problem,  optimal  control  problems  with 
solutions  published  in  the  literature  have  solutions  with  few  switches 
in  the  control.   The  multiple  switches  in  the  control  function  and 
oscillations  in  the  state  trajectory  result  in  a  number  of  local 
minima.   A  steepest  descent  algorithm  can  get  stuck  in  a  loop  at  a 
local  minimum.   For  the  Formulations  2  and  3  it  was  found  that  several 
local  minima  of  the  composite  cost  function  existed  which  did  not 
satisfy  the  optimality  conditions  of  the  original  problem.   In  such 
cases  it  was  necessary  to  exit  the  algorithm  and  make  manual  changes 
before  continuing  the  iterations. 
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Formulation  1  was  tried  initially  since  it  was  a  direct 
application  of  a  simple  gradient  method  to  the  solution  of  the 
original  problem.   The  basic  premise  of  the  controllability  of  the 
variational  system  naturally  arises  herein.   A  first-order  approximation 
of  an  increment  in  the  final  state  due  to  increments  §u   and  ou^  j.n 
the  control  and  dt   in  the  final  time  is  given  by  (cf .  ,  Equation  (3.8.11b)) 

dX(tf)  =  f(tf)dt  +  JQ  RTf ,u   5ux  dt  +  JQ  RTf,u  5u2  dt  . 

Cl      1  C2      2 

If   dt  =  0 


dX 


<V  =  L   ^'u     H  dt  +  L  ^-'u      6u2  dt  •  (3-9-1} 


Cl       1  C2 


It  is  assumed  that  5u   and  5u   are  piecewise  continuous. 

Let  the  control  increments  be  the  linear  combinations 

6u-  (t)  =  vT(t)  v   and   ftu  (t)  =  V  (t)   v  (3.9.2) 

1        —  X     —  ^        ~*^ 

where  V    (t)    and  V    (t)    are   suitable  vectors   of   dimension  6    and  y   is 

a   constant   vector   of   dimension  6.       Equation    (3.9.1)    may  be  rewritten   as 

dX(t    )    =  [  4   RTf,u      £    dt  +   Jo   RTf,u      f2   dt]  v    .  (3.9.3) 


It  may  be   seen  that    a   small    arbitrary   change  dX(tf)    in   the  final    state 
is   possible   if   and  only    if   the  matrix 

f      RTf,         Vjdt+J     RTf.         V^    dt]  (3.9.4) 

|_  Jo        -'u      -1  Jo         -   u      -2        J 


Cl  X  C2 


is    nonsingular. 
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Theorem 


There  exist  V  (t)  and  V  (t)  such  that 


r>        T       T       f>   T       T 

R  f ,    Vdt+     Rf,    Vdt 
«J  O    -'u   -         <Jq    -  U_  - 


T       T       0   T       T 
V 

1   1 


'1       _   ~        2 
is  nonsingular  if  and  only  it 

P  RTf,    fT,        R  dt  +  f   RTf,    fT    R  dt  (3.9.5) 

1  £ 

is  nonsingular. 

Proof   (Extension  of  Herme' s  [34]  theorem  (1)) 

T  T 

"if  part"  is  obvious  if  one  chooses  V  =  R  f,    and  V  =R  f,   . 


1        1 

To  prove  the  "only  if  part,"  assume  that  the  matrix 

n         T        T  •»    T       T 

R  f .    f,    R  dt  +  i   R  f,    f,    R  dt 
Jo    -  ux  ->Ul        Jco    -  u2  -  u2 

1  ^ 


-2 


(3.9.6) 


is  singular  but  the  expression  (3.9.4)  is  nonsingular. 

If  the  matrix  of  the  expression  (3.9.6)  is  singular,  then  there 
exists  a  vector   c  t-   0  such  that 


"" '  rp  rn  ft  T  T 

I   R  f,    f,   R  dt  +    R  f ,    f ,   R 
Jo   -  u_  -  u„        io    -  u„  -  u0 


LC! 


"i    "i 


T 
■'U2 


dt 


c  =  0    (3.9.7) 


or 


J' 


T      T 
R  f ,    f ,    R  dt 

"  Ul  -  Ul 


c  +  c 


T 
R  f, 


u„  -  u„ 


R  dt 


Since  both  the  above  terms  are  positive  semidef inite,  they  must 
be  individually  zero.   Also,  the  integrands  of  the  above  expressions 
are  positive  semidef inite  for  all  time.   So,  we  finally  have 
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T  T  o    "\ 

c  R   f ,        =0        almost    everywhere  on   C, 

"  Ul  " 


(3.9.8) 


T  T  ° 

c  R  f ,    =0    almost  everywhere  on  C 

-  u2  2  J 


Now,  premultiplying  the  expression  (3.9.4)  by  c   and  post- 
multiplying  by  c  and  using  the  results  (3.9.8),  one  obtains 


~        rp  rp  .       rp  rn  — 

R  f ,    V.,  dt  +     R  f ,    V   dt 
Jo    -  u   -1      <-'o    -  U„  -„ 


LC1 


2   2 


c  =  0 


or 


P  RTf,    V^  dt  +  f  R  f ,    V^  dt 
Jco   ->u±   -1      J,o   -  u2  -2 

is  singular. 

This  contradicts  the  assumption  that  this  matrix  is  nonsingular 
and  completes  the  proof  of  the  above  theorem. 

From  now  on  the  matrix  (3.9.6)  will  be  called  the  local  con- 
trollability matrix,  or,  simply,  the  controllability  matrix. 

In  all  solution  attempts  with  Formulation  1,  this  matrix  was 

-16      -19 
found  to  have  determinants  of  the  order  of  10    to  10    with  values  of 

the  elements  ranging  from  0.1  to  160.0.   The  computed  values  of  6u 
were  of  the  order  100  ft-lb  in  some  intervals.   The  large  §u  values 
persisted  even  when  small  changes  in  the  terminal  errors  were  programmed 
with  no  change  in  the  final  time  allowed.   The  large  §u  caused  insta- 
bility in  the  iterations.   This  indicated  that  the  terminal  state  com- 
ponents cannot  be  reached  simultaneously  with  control  increments  of 
this  method  for  the  trajectory  constructed  previously.   Rather  than 
attempt  to  reduce  the  terminal  state  components  simultaneously,  it 
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was  decided  to  form  a  norm  of  the  terminal  error  to  improve  the  local 
controllability  for  reducing  the  terminal  error.   This  approach  was 
used  in  Formulations  2  and  3. 

The  cost  functional  chosen  for  minimization,  the  final  time,  was 
another  source  of  difficulty.   At  the  minimum  time,  the  terminal  states 
are  just  barely  reachable.   The  exact  minimum  time  solution  cannot  be 
obtained  by  the  steepest  descent  algorithms  like  Formulation  1,  because 
the  controllability  matrix  becomes  singular  for  the  minimum  time  tra- 
jectory.  (Otherwise,  a  small  negative  dt  may  be  prescribed  for  which 
suitable  change  in  the  control  could  be  found  which  would  bring  the 
terminal  error  at  t  -  dt   to  zero.   This  would  contradict  the  minimality 
of  t  . )   Computational  difficulties  with  inversion  of  this  matrix  may  be 
expected  to  appear  some  time  before  the  actual  minimum  time  is  reached. 
The  main  problem  was  found  to  be  the  determination  of  the  increment  in 
the  final  time.   The  final  time  itself  being  the  cost  function,  it  is 
desired  that  the  increment  in  the  final  time  be  negative.   On  the  other 
hand,  convergence  of  the  terminal  errors  to  zero  may  require  an  increase 
in  the  final  time.   Unless  proper  weights  are  given  to  these  tendencies, 
the  value  of  the  final  time  may  become  either  larger  or  smaller  than 
the  actual,  minimum  time  of  the  problem.   Which  one  of  the  two  has 
occurred  cannot  be  determined  as  long  as  the  terminal  errors  are  not 
small.   If  the  errors  are  acceptably  small,  there  would  remain  the  pos- 
sibility that  the  final  time  can  be  reduced  further.   On  the  other  hand, 
if  the  final  time  has  been  unknowingly  reduced  to  less  than  the  minimum 
final  time  of  the  problem,  and  emphasis  is  greater  on  the  reduction  of 
the  final  time,  the  terminal  errors  will  remain  large.   Thus,  the 
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minimum-time  solution  for  a  complicated  problem  is  likely  to  result  in 
considerable  computations  with  these  methods. 

During  computations  with  Formulations  2  and  3,  it  was  found  that 
local  minima  of  the  composite  cost  functional  existed.   Rapid  improve- 
ment of  the  norm  of  the  terminal  error  was  made  in  the  first  few  iter- 
ations from  a  poor  control  function  guess;  but  the  rate  of  convergence 
became  extremely  slow.   The  changes  in  the  control  and  the  state  tra- 
jectories became  extremely  small  between  successive  iterations  once  the 
errors  became  small.   Formulations  2  and  3  do  not  guarantee  reduction  of 
all  the  terminal  errors  simultaneously  for  a  given  set  of  the  weights  K. 
and  it  was  therefore  necessary  to  readjust  these  weights  to  reduce  all 
the  errors  in  a  reasonable  number  of  iterations. 

The  important  events  during  the  computations  are  now  given. 
Formulation  2  was  used  first  due  to  the  following  reasons: 

a.  It  is  more  difficult  to  make  a  good  guess  for  a  constrained- 
unconstrained  arc-type  control  than  for  a  control  which  is  free 
to  change  everywhere. 

b.  In  Formulation  2,  the  changes  in  the  control  are  computed  every- 
where on  the  trajectory,  but  in  Formulation  3  the  changes  on 
the  unconstrained  parts  are  estimated  by  linear  approximations 
for  small  changes  in  the  state  variables.  It  was  found  that  the 
change  in  the  constraints  were  large  due  to  small  changes  in  the 
states.  Therefore,  the  algorithm  of  Formulation  3  was  likely  to 
be  less  stable,  especially  when  started  from  a  poor  guess. 

c.  Formulation  2,  after  it  ceases  to  give  improvement,  was  expected 
to  result  in  a  good  initial  guess  for  starting  Formulation  3. 
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Improvement  was  measured  by  the  trend  of  the  control  to  a 
bang-bang  control  or  a  stable  of f-the-constraint  arc  (for 
singular  subarc)  with  a  simultaneous  reduction  in  the  terminal 
error. 

Formulation  2  started  with  a  nominal  guess  shown  in  Figure  13, 
with  the  terminal  error  |  =  [1.66,  7.55,  -2.12,  -15.3,  -0.237,  22.5]  . 
It  may  be  noticed  that  the  final  state  in  the  initial  guess  is  far 
away  from  the  target  point.   The  model  has  reached  only  half  the 
required  height.   It  was  therefore  anticipated  that  the  nominal  final 
time  was  less  than  the  minimum  final  time  for  the  problem.   Thus,  by 
putting  less  emphasis  on  the  reduction  of  the  final  time  (keeping  Kq 
small)  the  minimum  time  solution  would  be  approached  monotonically  with 
the  decrease  in  terminal  error. 

A  local  minimum  of  the  cost  function  was  obtained  after  50  iter- 
ations.  This  was  determined  by  examination  of  the  gradient,  which  was 
quite  small.   The  quantity  B  (defined  in  Equation  (3.8.33)),  the  integral 

norm  of  the  gradient  of  the  cost  functional  with  respect  to  the  controls 

3 
had  a  value  of  0.954.   This  quantity  was  1.11  X  10   for  the  initial  guess. 

The  terminal  error  was  not  quite  acceptable  and  the  control  was  off  the 

boundary  most  of  the  time.   The  control  was  changed  manually  by  adding 

another  switch  in  the  control  near  the  final  time. 

After  75  more  iterations,  computations  were  stopped  when  the 

terminal  error  |  was  reduced  to  [-0.00124,  -0.0188,  -0.0834,  0.0993, 

T 
-0.0819,  0.391]   and  the  convergence  rate  was  very  poor.   The  control 

at  this  stage  is  shown  in  Figure  14.   The  final  time  was  2.0268  seconds, 

a  reduction  of  1.23  percent  from  the  initial  guess.   It  was  noticed 
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that  the  constraint  boundaries  had  moved  far  away  from  their  original 
locations.   On  the  other  hand,  all  the  significant  changes  in  the  con- 
trol had  been  small  compared  to  u  and  occurred  in  the  first  15  iterations. 

Since  the  optimal  solution  should  be  bang-bang,  or,  possibly, 
bang-bang  with  singular  subarcs,  the  control  was  modified  again  by 
placing  the  controls  on  the  nearer  of  the  two  boundaries.   This  was 
done  in  an  attempt  to  reduce  the  final  time  by  using  the  apparent  excess 
control  energy. 

These  changes  in  the  control  initially  increased  the  terminal 
errors  but  also  increased  the  prospects  of  meeting  the  terminal  condi- 
tions in  slightly  less  time.   The  terminal  error  decreased  in  the  sub- 
sequent iterations.   This  time  it  was  found  that  over  several  segments, 
especially  on  the  second  half  of  the  control  interval ,  the  constraint 
boundary  moved  inward  but  the  increments  in  the  control  were  computed 
outward.   As  a  result,  the  actual  change  in  the  control  used  when  the 
control  was  not  allowed  to  violate  the  constraints  increased  the  cost 
functional  instead  of  decreasing  it.   If  the  controls  were  allowed  to 
violate  the  constraints,  the  cost  decreased. 

At  this  point  it  was  noted  that  the  terminal  errors  were  quite 
sensitive  to  small  changes  in  the  control.   The  final  conditions  could 
be  met  by  adjusting  the  controls  only  in  the  later  parts  of  the  tra- 
jectory.  Slight  changes  in  the  earlier  parts  of  the  control  history 
caused  two  problems.   First,  they  caused  significant  changes  in  the 
control  boundaries  which  were  in  directions  opposite  to  the  directions 
for  improved  cost.   Second,  the  computed  control  increments  were  less 
accurate  near  the  beginning  of  the  solution  time,  so  that  with  high 
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sensitivity,  improvements  were  less  likely  via  control  increments  near 
the  beginning.   When  no  changes  were  allowed  in  the  earlier  part  of  the 
control  history,  the  terminal  errors  reduced  to  [0.0970,  -0.136,  -0.0883, 
-0.0340,  -0.0790,  -0.0438]  at  the  17  5   iteration.   The  final  time  was 
1.992  seconds.   The  convergence  rate  of  the  cost  function  was  again 
extremely  slow  in  the  last  several  iterations.   The  gradient  of  the  cost 
functional  with  respect  to  the  control  was  small  and  the  gradient  with 
respect  to  a   was  found  to  oscillate  between  positive  and  negative  values 
during  the  last  several  iterations.   The  control  history  is  shown  in 
Figure  15. 

This  control  was  considered  to  be  very  nearly  the  minimum  time 
control.   From  the  solution  it  was  clear  that  there  were  no  singular 
subarcs  and  the  optimal  control  would  be  bang-bang  if  further  refine- 
ment of  the  solution  was  made.   However,  observing  the  sensitivity  of 
the  final  time  to  the  control  changes  from  Figure  14  to  Figure  15  (about 
1%  for  large  control  changes)  it  was  believed  that  no  significant  reduc- 
tion in  the  final  time  could  be  found. 

Formulation  3  was  used  to  attempt  further  improvements  in  the 
solution  and  to  check  the  validity  of  the  above  assertions.   The  con- 
trols were  modified  by  placing  them  exactly  on  the  constraints  when 
they  were  close  to  the  constraints,  and  off  the  constraint  arcs  were 
made  steeper.   The  terminal  errors  initially  increased  because  of  these 
changes.   Formulation  3  was  found  to  converge  more  slowly  than  Formula- 
tion 2.   This  was  due  to  the  fact  that  in  Formulation  3  the  changes  in 
the  control  were  made  on  relatively  small  segments.   The  changes  in  the 
constraints  were  found  to  be  larger  than  those  of  the  control  increments 
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in  the  unconstrained  parts.   This  was  because  the  constraints  were  very 
sensitive  to  changes  in  the  states. 

The  algorithm  for  Formulation  3  was  able  to  reduce  the  terminal 
errors  for  90  iterations  until  it  was  necessary  to  stop  increments  in 

the  first  half  of  the  control  function.   Again  the  gradient  there  became 

th 
noisy  and  inaccurate.   Computations  were  stopped  at  the  120    iteration 

when  the  errors  were  [0.0504,  0.869,  -0.0366,  0.0154,  0.0081,  -0.113], 
The  errors  for  the  state  variables  X^,  Xg,  X^,    Xg,  and  XQ   were  accept- 
ably small  and  the  error  for  Xg  was  slightly  greater  than  what  was 
obtained  before  by  Formulation  2.   The  values  of  the  6u  in  the  final 
ten  iterations  were  of  the  order  0.01  ft-lb.   Larger  values  resulted 
in  divergence  or  an  increase  in  cost. 

The  control  history  obtained  by  Formulation  3  is  shown  in 
Figure  16.   The  final  time  in  the  trajectory  was  1.992,  the  same  as 
the  final  time  of  Figure  15.   The  final  time  was  less  in  the  intermedi- 
ate iterations,  but  it  increased  back  to  this  value.   From  these  results 
it  was  concluded  that  the  solution  should  indeed  be  bang -bang,  but  no 
significant  improvement  in  the  final  time  could  be  achieved. 

It  may  be  pointed  out  that  the  solution  (b)  of  the  Example 
Problem  in  Section  3.7,  shows  the  final  time  2.025  was  only  1.25  percent 
more  than  the  exact  bang-bang  solution  final  time  of  2.000.   But  in  the 
solution  (b)  the  control  was  off  the  constraint  (during  the  switch)  for 
28  percent  of  the  total  time. 

A  preliminary  investigation  of  a  steepest  descent  scheme  for 
bang-bang  solution  is  described  in  Appendix  B. 
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The  difference  between  the  state  trajectories  for  the  solutions 
of  Figures  15  and  16  was  insignificant.   The  trajectories  for  the  angles 
cp,  6,  and  J  are  shown  in  Figure  17  for  the  solution  of  Figure  15. 

Examination  of  the  Computed  Solution 
by  the  Principles  of  Mechanics 

The  computed  solution  can  be  explained  qualitatively  by  consider- 
ing the  average  motion  of  the  links  and  the  mechanics  of  rotating  systems. 
We  notice  that  positive  controls  tend  to  reduce  the  radius  of  gyration 
of  the  system  at  most  configurations  and  hence  they  tend  to  increase 
jcpj  which  is  the  major  contributor  to  the  angular  momentum  of  the 
system.   The  negative  controls  do  the  opposite. 

The  controls  remain  positive  during  most  of  the  time  in  order 
to  obtain  increased  jcpj  except  near  the  end  of  a  swing,  i.e.,  at  the 
maximum  amplitude.   There,  the  controls  switch  to  negative  values  in 
order  to  limit  the  swing  prior  to  the  final  swing  and  thus  save  time. 
It  should  be  noted  that  the  first  swing  should  be  just  high  enough  so 
that  with  his  limited  strength,  the  gymnast  is  able  to  swing  to  the  top 
on  the  second  swing.   Thus,  during  the  first  two  times,  the  angle  9 
reached  its  maximum  amplitudes.   There  are  four  switches  each  in  the 
controls  u  and  u  .   At  the  end  of  the  last  swing  we  should  not  expect 
any  more  switch.   However,  in  order  to  meet  the  terminal  conditions, 
the  legs  had  to  be  brought  down,  requiring  another  switch  in  u2- 
This  switch  also  helped  the  angles  9  and  cp  to  reach  their  terminal 
values.   After  this  has  happened,  one  more  switch  each  for  the  controls 
was  needed  to  meet  the  terminal  conditions  on  the  angular  rates. 
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Without  any  terminal  conditions  on  the  angular  rates,  the  solution 
could  have  been  obtained  without  the  last  switch  in  each  of  the  controls 
and  in  slightly  less  final  time. 

3.10.   Comparison  of  the  Minimum-Time 
Solution  With  Experiment 

A  comparison  between  the  control  histories  of  the  computed 
solution  and  the  maneuver  of  the  gymnast  recorded  in  the  experiment 
with  the  same  boundary  conditions  was  not  possible  because  the  control 
forces  used  by  the  gymnast  could  not  be  determined.   The  photographic 
measurements  of  the  angles  were  not  accurate  enough  for  two  successive 
numerical  differentiations  needed  for  finding  the  control  torques  for 
the  actual  rigid  body  motion  from  the  equations  of  motion.   The  time 
histories  of  the  generalized  coordinates  cp ,  9,  and  f  for  the  experiment 
and  the  computed  solution  are  given  in  Figures  5  and  17.   The  trajec- 
tories of  9  and  l|f  are  quite  different  in  the  two  cases.   The  differ- 
ences are  primarily  due  to  the  differences  between  portions  of  the  con- 
trol limit  functions  used  in  the  mathematical  model  and  the  actual 
strength  of  the  gymnast.   During  the  simulation  runs  it  was  found  that 
the  model's  arm  strength  was  insufficient  for  small  values  of  8.   Thus, 
the  actual  state  trajectory  was  not  a  good  guess  for  starting.   A  nom- 
inal trajectory  had  to  be  determined  which  was  more  nearly  correct  for 
a  gymnast  with  less  strength  at  extremes  of  his  limb  movement. 

It  may  be  noticed  that  in  the  measurements  of  the  angle  9  of 
the  kip-up  motion,  the  angle  9  remained  non-negative  throughout  the 
motion,  whereas  in  the  computed  minimum  time  solution  8  became  negative 
for  a  short  duration,  with  a  minimum  value  of  -.25  radians.   This 
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discrepancy  is  largely  due  to  inaccuracy  in  modeling.   The  modeling 
inaccuracies  include  modeling  the  head-torso  system  as  a  rigid  link, 
and  inaccuracies  in  the  inertia  parameters  and  the  functions  describ- 
ing the  constraints  on  the  control  torques. 

From  the  physical  considerations,  it  can  be  seen  that  raised 
arms  can  indeed  bend  more  than  .  25  radians  behind  the  torso.   But  then 
the  bending  is  shared  by  the  deformation  of  the  torso,  as  demonstrated 
in  Figure  18.   Under  such  conditions  there  will  be  differences  in  the 
angles  8  and  i|f  measured  from  the  tape  strips  and  the  corresponding 
angles  of  the  model  based  on  lines  between  joint  centers.   This  is 
illustrated  in  Figure  18. 

In  the  optimization  procedure,  an  inequality  constraint  on  the 
angle  9  would  be  able  to  keep  8  from  being  less  than  a  prespecified 
value.   Such  a  constraint  was  not  considered  in  this  investigation  for 
several  reasons.   First,  it  was  thought  that  if  the  strength  measure- 
ments (the  control  limits)  were  accurate,  the  states  would  remain  within 
acceptable  bounds  naturally  and  without  requiring  a  constraint  on  the 
states  imposed  externally.   Second,  the  exact  limits  and  the  nature 
(hard  or  soft)  of  the  constraints  on  the  angles  need  further  exper- 
imental investigations.   And,  lastly,  it  was  anticipated  that  the 
violation  of  the  state  constraint,  if  any,  would  not  be  significant 
enough  to  justify  going  into  the  complications  associated  with  hard 
state  constraints. 
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Figure  18.   Difference  Between  Measured  Angle  and  Mathematical 
Angle  for  Human  Model  Due  to  Deformation  of  Torso. 


CHAPTER  4 


CONCLUSIONS  AND  RECOMMENDATIONS 
FOR  FUTURE  WORK 


4.1.   Conclusions 

The  results  of  an  investigation  involving  application  of 
rational  mechanics  and  optimal  control  theory  to  a  human  motion  prob- 
lem have  been  presented.   The  kip-up  maneuver  of  gymnastics  was 
selected  for  this  purpose.   First,  a  three-link  model  was  constructed 
for  the  human  performer — a  professional  gymnast.   The  dynamic  prop- 
erties of  the  model  were  compared  with  those  of  the  gymnast  by  means 
of  experiments.   A  minimum-time  strategy  for  the  kip-up  maneuver  was 
then  obtained  by  numerical  computation  for  the  mathematical  model. 
The  solution  obtained  was  compared  with  an  actual  minimum  time  maneu- 
ver performed  by  the  gymnast. 

It  was  found  that  the  human  being  may  be  modeled  by  rigid  links 
for  the  purpose  of  dynamic  analysis  of  human  motion  and  man-machine 
systems.  However,  the  model  cannot  be  tested  effectively  unless  accur- 
ate measurements  of  the  angular  velocities  and  angular  accelerations 
along  with  the  angles  are  made.  Photographic  determination  of  the  limb 
positions  is  not  accurate  enough  for  determining  within  the  uncertain- 
ties of  the  model  the  derivatives  of  the  movements  involved. 
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Modifications  were  needed  for  the  shoulder  and  hip  joints  to 
account  for  the  stiffness  and  the  deformation  of  these  joints  during 
extreme  arm  and  leg  movements. 

An  accurate  knowledge  of  the  maximum  strength  of  the  person 
being  modeled  at  the  junctions  of  the  various  body  segments  is  neces- 
sary for  correctly  determining  the  optimal  performance.   For  this 
reason,  these  measurements  must  be  taken  at  small  intervals  of  the 
joint  angles  and  for  the  entire  range  of  the  possible  movements  so  that 
a  polynomial  representation  or  interpolation  of  these  limits  with 
respect  to  the  angles  is  sufficiently  accurate.   In  the  present  inves- 
tigation, the  control  limit  functions  used  did  not  represent  the  true 
strength  of  the  gymnast  at  the  extremities  of  the  arm  movement,  and  as 
a  result,  a  significant  difference  in  the  actual  maneuver  of  the  gymnast 
and  the  optimal  motion  of  the  mathematical  model  was  obtained. 

After  an  accurate  mathematical  model  has  been  constructed,  the 
optimal  motion  has  to  be  determined  by  numerical  computation. 

The  analytical  determination  of  the  minimum  time  kip-up  maneuver 
involved  working  with  a  significantly  nonlinear  system  with  control 
constraints  which  depend  on  the  state  variables.   Several  numerical 
methods  were  used  to  obtain  the  solution.   A  modified  quasilineariza- 
tion  method  was  tried  first  but  the  attempt  was  unsuccessful.   First- 
order  steepest  descent  methods  were  then  used  with  three  different 
formulations.   A  nearly  optimal  numerical  solution  was  obtained  from 
these  formulations.   The  solution  was  compared  with  experimental  results. 
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It  was  found  that  the  method  of  quasilinearization  was  not 
suitable  for  the  highly  nonlinear  human  motion  problems  involving 
multiple  switches  in  the  control.   It  was  also  noted  that  the  method 
would  not  be  efficient  for  problems  which  may  admit  singular  solutions. 
It  (or,  any  other  method  which  determines  the  optimal  solution  via 
solution  of  a  two-point  boundary  value  problem  of  the  state  and  adjoint 
variables)  uses  the  optimal  expression  for  the  control  variables  in 
terms  of  the  state  and  adjoint  variables.   Several  problems  would  be 
faced  in  determining  the  optimal  control  when  singular  arcs  are  to  be 
considered.   First,  the  equalities  in  the  necessary  conditions  for 
singular  subarcs  would  be  difficult  to  verify  in  the  limited  precision 
of  digital  computations.   Second,  sufficiency  conditions  for  optimality 
of  singular  arcs  for  a  general  nonlinear  system  have  not  yet  been  estab- 
lished.  Lastly,  if  the  initial  guess  of  the  state  and  adjoint  variables 
is  poor,  then  it  may  happen  that  in  an  intermediate  iteration  the  state 
and  adjoint  variable  histories  correspond  to  a  singular  control  arc  on 
a  certain  time  interval,  but  they  do  not  do  so  for  the  same  time  interval 
in  the  next  iteration.   In  such  a  case  the  iterations  are  likely  to  be 
unable  to  improve  the  solutions. 

The  method  of  quasilinearization,  or,  likewise,  any  second-order 
method  requires  more  computation  and  storage  than  first-order  methods. 
For  the  human  motion,  the  method  of  quasilinearization  required  an  exces- 
sive amount  of  computing  time  per  iteration  and,  more  critically,  nearly 
maximum  storage.   The  second-order  methods  were  therefore  not  considered 
further. 
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Existing  first-order  steepest  descent  methods  are  not  efficient 
for  solving  such  problems  either.   Though  the  solution  was  obtained  by 
these  methods,  simulation  and  physical  reasoning  played  an  important 
role.   A  Davidon-Fletcher-Powell  scheme  of  a  modified  gradient  method 
was  also  used  but  without  improvement  in  the  results. 

The  method  of  constrained  and  unconstrained  arcs  is  not  effec- 
tive as  a  starting  method  when  the  control  constraints  of  the  constrained 
and  unconstrained  arcs  are  not  known  in  advance.   The  effectiveness  of 
the  constrained-unconstrained  arc  method  depends  on  the  sensitivity  of 
the  control  constraints  to  changes  in  the  state  variables.   Its  conver- 
gence is  poor  when  the  constraints  are  sensitive  to  changes  in  the 
states.   On  the  other  hand,  treating  the  control  as  entirely  free  while 
computing  the  gradients  and  later  imposing  the  constraints  while  imple- 
menting the  control  is  a  much  simpler  scheme  and  converges  more  rapidly 
when  a  significant  portion  of  the  control  history  is  away  from  the  con- 
straints.  The  method  becomes  ineffective  (and  theoretically  incorrect) 
when  controls  are  mostly  near  the  constraints  and  the  movement  of  the 
constraints  becomes  significant. 

Each  of  the  steepest  descent  algorithms  requires  a  good  initial 
guess  when  several  switches  occur  in  the  optimal  control.   They  may  not 
be  able  to  produce  a  required  switch  in  the  control  in  problems  with 
multiple  control  switches. 

In  Formulations  2  and  3  of  the  steepest  descent  method,  the 
adjustment  of  the  final  tine  was  done  by  adjusting  the  state  variable  a 
introduced  by  changing  the  independent  variable  from  t  to  t ,  where  t=  qt. 
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This  change  was  expected  to  improve  the  performance  of  the  steepest 
descent  method  with  free  final  time  because  a   appeared  in  the  state 
equations.   This  approach  was  not  found  to  be  more  efficient  than  the 
conventional  method  of  changing  the  final  time,  i.e.,  by  extending  or 
truncating  the  final  end  of  the  trajectory.   The  gradient  of  the  norm 
of  the  terminal  error  with  respect  to  a  was  extremely  sensitive  with 
respect  to  changes  in  a   as  well  as  other  changes  in  the  trajectory. 
This  required  the  change  6a  be  kept  extremely  small  in  iterations, 
using  nearly  optimal  results.   It  was  significant  that,  during  the 
first  manual  change  of  the  control  after  achieving  a  local  minimum, 
it  was  necessary  to  increase  a   to  its  original  value  of  1.0  (it  had 
dropped  to  0.977)  and  truncate  the  trajectory.   This  adjustment  is  the 

same  as  the  conventional  method  of  adjustment  of  the  final  time. 

5 
For  the  relatively  quick  kip-up  maneuver,  1.36  x  10  bytes 

of  storage  (with  double  precision  arithmetic)  and  a  computing  time  of 

about  six  seconds  per  iteration  (for  IBM  370-165  Computer,  Fortran  G) 

were  required  for  Formulations  2  and  3  of  the  steepest  descent  method. 

4.2.   Recommendations  for  Future  Work 

The  following  areas  need  to  be  considered  for  future  work  in 
the  study  of  dynamics  and  analytic  determination  of  human  motions: 
1.   Improvement  of  the  Inertia  Model  of  Hanavan  and  Data 
Gathering  Experiments: 

The  most  significant  inaccuracy  in  the  three-link  inertia 
model  for  the  kip-up  maneuver  was  observed  to  be  the  deformation  of  the 
shoulder  and  bending  of  the  torso  during  certain  periods  of  the  motion. 
An  accurate  determination  of  the  average  location  of  the  joint  centers 
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at  the  shoulder  and  hip  is  sufficient  for  modeling  the  middle  link  of 
the  three-link  kip-up  model. 

The  stiffnesses  of  the  hip  and  shoulder  joints  during  extreme 
arm  and  leg  movements  should  also  be  determined  by  experiment. 

The  bending  of  the  torso  might  be  accurately  modeled  by 
dividing  the  torso-head  link  into  two  links.   The  upper  link  might 
consist  of  the  head  and  upper  torso  with  the  lower  torso  of  the 
Hanavan  model.   These  two  links  might  be  joined  by  a  smooth  hinge  with 
a  torsional  spring  and  a  voluntary  torque.   The  introduction  of  the 
additional  degree  of  freedom  to  take  care  of  the  deformation  of  the 
torso  will  increase  the  complexity  of  the  dynamical  equations  of  the 
system  further. 

Work  needs  to  be  done  to  improve  the  data  gathering  techniques, 
that  is,  determination  of  the  generalized  coordinates  and  their  rates. 
Photographic  measurements  are  extremely  time  consuming  and  inaccurate. 
Preliminary  evaluations  [35]  show  that  the  angular  differentiating 
accelerometers  (ADA)  will  provide  the  desired  rate  and  angle  accur- 
acies.  The  devices  are  quite  small  and  have  a  low  power  requirement. 
With  a  light  weight  telemetry  system  the  human  performer  is  not  encum- 
bered by  wiring.   Once  accurate  determination  of  the  angles  and  their 
rates  is  possible,  the  inertia  model  may  be  improved  by  using  the 
experimental  data  and  some  of  the  system  identification  techniques. 
2.   Measurement   of  the  Muscle  Torques: 

An  accurate  knowledge  of  the  strength  of  the  individual  is 
needed  in  the  analytic  determination  of  a  given  maneuver.   It  needs  to 
be  investigated  how  the  strength  at  a  joint  depends  on  the  rate  of 
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contraction  of  the  muscle  group  and  the  time  for  which  the  force  acts. 
The  actual  muscle  forces  or  their  equivalent  torques  during  a  maneuver 
should  also  be  determined  in  order  to  have  a  better  understanding  of 
the  muscle  activities.   If  accurate  information  of  the  inertia  proper- 
ties and  the  angles  and  their  rates  is  available,  these  forces  can  be 
obtained  from  the  equations  of  motion.   However,  determination  of  these 
forces  using  biomedical  devices  such  as  electromyography  electrodes 
might  give  valuable  information  in  this  area. 

3.   Determination  of  Better  Numerical  Techniques  for  Human 

Motion  Problems: 

It  was  clear  from  the  present  investigation  of  the  two- 
second  duration  kip-up  maneuver  that  analytic  solution  of  more  compli- 
cated and  longer  duration  motion  would  be  extremely  difficult.   Consider- 
ation of  state-constraints  may  be  necessary  to  ensure  that  the  stresses 
generated  during  the  motion  at  any  part  of  the  human  body  remain  within 
certain  limits.   Also,  constraints  should  be  imposed  on  the  rate  of 
change  of  the  control  forces  because  a  human  being  may  not  be  able  to 
execute  bang-bang  control.   These  constraints  will  further  complicate 
the  optimization  problems  and  introduce  interesting  biomechanical 
questions. 

It  therefore  seems  that  some  numerical  methods  should  be  found 
which  would  give  an  approximate  extremum  (i.e.,  not  give  an  exact 
extremum)  but  would  be  efficient  in  handling  the  constraints  and  other 
complications  of  the  human  motion  problem. 
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APPENDIX  A 


DETERMINATION  OF  THE  INERTIA  PARAMETERS 
OF  THE  KIP-UP  MODEL  FROM  THE  HANAVAN  MODEL 


The  three  elements  of  the  kip-up  model  were  constructed  by 
combining  the  elements  of  the  Hanavan  model.   The  inertia  properties 
of  the  kip-up  model  were  obtained  from  geometric  considerations  and 
the  parallel  axes  and  rotation  theorems  of  moment  of  inertia. 
Figure  19  shows  the  three  elements  of  the  kip-up  model  and  their  com- 
ponents from  the  Hanavan  model.   The  dimensions  of  the  Hanavan  model, 
including  the  location  of  the  hinge  axes  A  and  B,  could  be  obtained 
from  the  program  of  Hanavan.   Input  of  this  program  were  the  25  anthro- 
pometric dimensions  taken  from  the  subject.   The  formulas  for  the  var- 
ious dimensions  of  the  elements  of  the  kip-up  model  are  presented  below, 
These  formulas  were  added  to  Hanavan' s  program  in  order  to  obtain  the 
inertia  parameters  of  the  kip-up  model. 

SM(8)[R(4)+fl-T](8)}SL(8)]cos  y +  SM(6)  [R(4)  +  SL(8)+f  1-Tj  (8)  }SL(6)  ]cos  y 
rl=  SM(4)  +  SM(8)  +  SM(6) 

I     =  [R(4)  +  SL(8)  +  SL(6)  -  R(6)]  cos  y 

m  =  2[SM(4)  +  SM(8)  +  SM(6)] 

-I   =  SIYY(4)  +  SM(4)rJ+^[SIYY(8)  +  SIZZ(8)  +  (SIYY(8)  -SIZZ(8)} 

cos  2y]  +  SM(8)[r1-[R(4)  +  SL(8)  -  71(8)  SL(8)  }cos  y]  +-[SIYY(6) 
+  SIZZ(6)  +  fSIYY(6)  -SIZZ(6)]cos  y]  +  SM(6)  [f  R(4)  +  SL(8) 


+  SL(8)  +  SL(6)  -  ~(6)}cos  y  -  r^\ 
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SLCiO) 


SLC2) 


Figure  19.      Construction  of  Kip-Up  Model   from  Hanavan  Model. 
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-SM(1)R(1)  +  SM(2)T)(2)SL(2)  +  SM(3)  [SL(2)+T]  (3)  SL(3)  ] 
P2  ~  SM(1)    +    SM(2)    SM(3)  (    } 

I     =    SL(2)    +    SL(3)    -   DELSH    -  R(6) 


m     =   SM(1)   +    SM(2)   +    SM(3) 


I2  =    SM(1)[R(1)   +   R(6)   +    r2]      +    SIYY(l)   +    SM(2)[R(6)   +    rg  -  Tl(2) 
SL(2)]2  +    SIYY(2)   +    SM(3)[R(6)   +   r     -   SL(2)    -  7]  (3)SL(3)]2 
+    SIYY(3) 
jG3  =   DELSH  +   |sM(10)T](10)SL(10)    +    SM(12)  [SL(10)    +    T)(12)SL(12) 

+    SM(14)[SL(10)    +    SL(12)+R(14)]j/|SM(10)    +    SM(12) 
+    SM(14)| 

m     =    2[SM(10)    +    SM(12)    +    SM(14)] 

•3 

i  Ip  =    SIYY(IO)    +    SM(10)[£3    -   DELSH    -  T](10)SL(10)]2 

.2 


+    SIYY(12)   +    SM(12)[i3    -   DELSH    -    SL(10)    -  T[  (12)SL(12)] 

*3 


+    SIYY(14)    +    SM(14)[SL(10)    +    SL(12)   +    R(14)    -   i     +    DELSH] 
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Anthropometric  Dimensions  of  the  Subject 

Dimension  Value 


(inches) 

1.  Ankle  Circumference  9.25 

2.  Auxiliary  Arm  Circumference  13.0 

3.  Buttock  Depth  8.125 

4.  Chest  Breadth  12.75 

5.  Chest  Depth  9.25 

6.  Elbow  Circumference  11.0 

7.  Fist  Circumference  11.0 

8.  Forearm  Length  (Lower  arm  length)  10.75 

9.  Foot  Length  10.0 

10.  Knee  Circumference  13.5 

11.  Head  Circumference  22.5 

12.  Hip  Breadth  12.0 

13.  Shoulder  Height  (Acromial  Height)  54.0 

14.  Sitting  Height  34.75 

15.  Sphyrion  Height  3.5 

16.  Stature  66.5 

17.  Substernal   Height  46.5 

18.  Thigh  Circumference  20.625 

19.  Tibale  Height  16.125 

20.  Trochanteric  Height  34.25 

21.  Upper  Arm  Length  12.0 

22.  Weight  138  lb 

23.  Waist  Breadth  10.25 

24.  Waist  Depth  7.625 

25.  Waist  Circumference  6.75 
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Inertia  Properties  of  the  Kip-Up  Model 


I     =   22.29  inches 
I     =      17.68  inches 


lengths  of  the  elements  1  and  2 


r  =  12.68  inches 


2.57  inches 
12.32  inches 


distances  of  the  CG' s  of  elements  1, 
2 ,  and  3  from  their  upper  hinge  points 


ml  = 

0.47 

slug 

m2  = 

2.39 

slug      / 

m3  = 

1.42 

slug    j 

I.  = 

28.20 

•  2^ 

slug-m 

208.98  slug-in 


140.90  slug-in  J 


mass  of  the  elements  1,  2,  3 


moments  of  inertia  about  the  centers 
of  gravity  for  elements  1,  2,  3 


APPENDIX  B 


AN  INVESTIGATION  OF  A  STEEPEST  DESCENT  SCHEME  FOR  FINDING 
OPTIMAL  BANG-BANG  CONTROL  SOLUTION  FOR  THE  KIP-UP  PROBLEM 


An  investigation  was  made  to  determine  the  suitability  of  the 
steepest  descent  algorithm  of  Bryson  and  Denham  [28]  for  determining 
optimal  bang-bang  control  for  the  kip-up  problem.   The  investigation, 
described  in  this  section,  showed  that  extremely  small  changes  in  the 
switching  times  were  necessary  for  the  first-order  relationship  between 
the  changes  in  the  switching  time  and  changes  in  the  terminal  state 
vector  to  be  valid.   The  small  changes  in  the  switching  time  ruled 
out  the  use  of  fixed  integration  step  size  as  used  for  the  continuous 
control.   The  change  of  the  switching  time  by  one  integration  step  was 
found  to  be  much  too  large  for  the  first-order  relationship  to  be  valid. 
Therefore,  an  attempt  to  find  a  bang-bang  optimal  control  by  using 
a  fixed  step  size  of  0.0036  second  (with  a  total  of  about  550  inte- 
gration steps)  was  not  successful.   Even  though  a  scheme  for  avoiding 
use  of  fixed  integration  step  size  was  available,  the  method  was  not 
pursued  further  because  of  anticipated  numerical  problems.   The  impor- 
tant results  of  the  investigation  are  now  presented  in  more  detail. 

Suppose  each  of  the  controls  u  (t)  and  u  (t)  is  on  either  the 

upper  limit  or  lower  limit  for  all  time.   Also,  let  the  control  u  (t) , 

i  =  1  or  2,    have  the  switching  times  t    ,t    , . . . ,t   ,  i.e.,  at  these  times 

i 
the  control  switches  from  an  upper  to  a  lower  limit  or  vice  versa. 
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We  shall  first  find  out  the  change  in  the  final  state  vector  due  to 
a  small  change  in  one  of  the  switching  times  t  .   For  this  we  can  use 
the  results  of  Section  3.8,  Equation  (3.8.11)  which  is  written  as 

6x(t  )  =  f(t  )dt+RT(0)6x(0)  +  J   RTf,    6udt  +  ifRTf,    au2dt-     t6-1) 
In  this  equation  C  and  C  denote  the  unconstrained  portions  of 

J.  Ct 

the  control  trajectories  a.  (t)  and  u  (t) ,  respectively.   In  the  present 
case  the  unconstrained  arcs  are  the  portions  of  the  control  trajectories 
which  change  from  one  constraint  to  another  due  to  the  shift  of  the 
switching  times.   They  are  therefore  around  the  switching  times  and  are 

not  known  a  priori.   But,  since  the  shifts  in  the  switching  times  will 

o 

be  made  small  in  any  iteration,  the  presence  of  the  integrals  over  C. 

and  C  can  be  ignored  as  second-order  terms  while  generating  R(t)  by 
backward  integration  of  Equations  (3.8.10).   The  value  of  R(t)  will 
therefore  be  the  solution  of  the  following  equations: 

R(t  )  =  1        (6x6  unit  matrix)  (B.  2) 

R  =  -{f   +  f.    S^    +  f     S*  f   R  (B.3) 

-    "  Ul   1'X    "  U2 


2'x 


where 


(B.4) 


j  =  1  for  u  on  the  lower  constraint 
j  =  2  for  iL  on  the  upper  constraint 
K  =  1  for  u  on  the  lower  constraint 
K  =  2  for  u   on  the  upper  constraint  .  _, 
Suppose  a  small  change  6x(0)  at  the  initial  condition  of  the 
state  vector  X(0)  is  prescribed  so  that  the  initial  state  becomes 
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X  +  §x(0).  Also,  let  the  switching  times  remain  unchanged,  and  the 
-o    - 

final  time  changes  from  t  to  t  +  dt  .  The  resulting  change  in  the 
final  values  of  the  state  variables  can  be  obtained  by  setting  &u =  0 
in  Equation  (B.l),  which  yields 


6x(tf)  =  R  (0)  6x(0)  +  f(tf)  dtf 


(B.5) 


If  at  any  other  time  t  ,  0  <  t  £  t  ,  the  stpte  vector  X(t  )  is 
changed  by  the  amount  §X(t  ),  the  change  in  the  final  state  vector, 
if  the  same  switching  times  are  maintained,  will  similarly  be  given  by 


6x(tf)  =  R  1(t1)  SXC^)  +  f(tf)  dtf 


(B.6) 


Now   consider  a  change  in  the  switching  time  t  .   Suppose  the 

2     1 
control  u   switches  from  S   to  S   at  this  switching  time,  and  that  this 

switching  time  will  be  advanced  by  the  amount  dt  .   All  other  switching 

times  remain  unchanged. 

Consider  the  value  of  the  state  vector  X(t)  at  the  time 

t  +  dt  .   Before  the  change  in  the  switching  time,  we  have,  to  the 

first  order 


-(tl+6tl)  =  -(tl}  + 


Ut\+) 


dt    (t   means  just  after  t  ) 


=  Xd^)  + 


x(tx)  + 


A  +  £'n   Ul  +  ^'uo  U2 


A  +  f ,    S^  +  f ,    u 
-  Ul  1         ~   U2   2 


u    < 


1  dtl 


(B.7) 


where  the  vector  A(X)  is  defined  in  Equation  (2.2.20). 
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After  the  change  in  the  switching  time,  we  have,  similarly, 


X(t*  +  6t*)  =  X(t±)   + 


A+  !>u   Sl  +  ^'u  U2 


(B.8) 


From  (B.7)  and  (B.8),  we  obtain  the  change  in  the  state  vector  at  time 


t  +  dt   due  to  the  change  in  the  switching  time. 


fix 


X(t*  +  6tJ)  =  f,u  (s2  -  sj)  dt* 


(B.9) 


12  12 

If  the  switch  was  from  S   to  S  we  would  have  obtained  (S1  -  S^) 

2    1  1 

instead  of  (S   -  S  )  in  Equation  (B.9).   If  d^  was  negative  the  rela- 
tion (B.9)  would  still  be  valid. 

Corresponding  to  the  change  in  the  state  vector  at  t  +  dt_^ 
given  by  Equation  (B.9),  the  change  in  the  final  state  vector  given  by 
Equation  (B. 6)  is 

6X(tf)  =  RT<ti+6ti>!'u  (S^-sJ)dtJ  +  f(tf)  dtf 


or,  to  the  first  order 


6X(tf)  =  f (tf)dtf  + 


rT£»u1<sX> 


idti 


(B.10) 


In  the  nominal  trajectory,  suppose  switches  from  S.  to  S. 

occur  at  t„  ,  t_,  tt , .  .  .    and  switches  from  S.  to  S.  occur  at 
13    5  ii 

t  ,  t  ,  t1 for  i  =  1  and  2.    The  result  (B.10)  can  be  extended 

2'   4'   5 

for  the  case  when  any  of  the  switching  times  t  is  advanced  by  dt  . 

This  will  yield 

N1 

i-1 


6x(t  )  =  f  (t  )dt  +  Z  (-D' 

i=l 


N2 
+   Z   (-D 
3=1 


j-l 


RTf,   (sj-sj) 
1 


*T*'u  (S2"S2) 
2 


1 
1  dti 


dt2 
.2    J 


(B.ll) 
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The  relations  (B.10)  and  (B.ll)  are  exactly  the  same  as  the 
gradient  relations  given  in  [24]  for  bang-bang  control,  but  are 
derived  differently.   This  derivation  brings  out  the  fact  that  the 
change  in  the  final  state  vector  will  be  given  by  the  linear  relations 
only  if  the  change  in  the  state  vector  at  a  switching  time  caused  by 
the  change  in  switching  time  is  small  enough  to  make  the  relation  (B. 6) 
valid. 

To  find  the  maximum  value  of  | 6X(t^) |  valid  for  use  in  expres- 
sion (B.6),  the  following  test  was  performed. 

First,  with  one  forward  and  one  backward  integration,  with 
a  nominal  trajectory,  the  value  of  R(0)  was  determined.   The  initial 
time  t  =  0  was  used  for  t  .   Several  values  of  SX(O)  were  chosen  and 
for  each  initial  state  vector  given  by  X  +  SX(0) ,  the  state  equations 
were  integrated  forward  with  the  same  switching  times  and  with  dt  =  0. 
The  final  values  of  the  state  vector  X(t  )  obtained  by  the  integration 
were  used  to  find  the  actual  changes  in  the  final  values  of  the  state 

vector  from  the  value  of  the  nominal  trajectory.   These  actual  changes 

T 
were  then  compared  with  the  predicted  change  R(0)  0X(0) . 

From  these  comparisons  it  was  found  that  6x(0)  with  elements  of 

the  order  of  0.01  was  too  large.   The  relation  (B.6)  was  valid  for 

a  much  smaller  (sy  0.001)  6x(0).   On  the  other  hand,  the  elements  of 

§X(t)  caused  by  a  change  in  the  switching  time  by  one  integration  step 

of  0.0036  second  used  in  Sections  3.7  and  3.8  were  found  to  be  of  the 

order  of  0.1  to  6.0 — much  larger  than  what  could  be  used  in  Equation 

(B.6).   From  the  matrix  R(0) ,  which  had  large  entries,  it  was  clear 

that  the  change  in  the  switching  time  should  be  much  smaller. 
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In  order  to  find  an  optimal  bang -bang  control  by  this  method, 
the  following  features  in  the  determination  of  the  influence  function 
R(t)  should  be  incorporated. 

1.  A  nominal  guess  consisting  of  the  switching  times  and  a  final 
time  is  to  be  made  first.   The  switching  times  specify  the 
control  completely  by  telling  which  constraints  the  control  is  on. 

2.  Integrate  the  state  equations  forward  with  the  initial  value  X  . 
Generate  u  from  the  constraints.   Store  X(t_).   It  is  not 
necessary  to  store  X(t)  and  u(t)  for  any  other  t. 

3.  Integrate  the  state  equations  backwards  and  the  influence 
Equations  (B.4)  simultaneously.   Use  final  value  of  the  state 
variables  as  those  obtained  from  the  forward  integration. 
Generate  control  from  the  constraints. 

Generating  the  state  and  the  control  variables  again  during 
the  backward  integration  (they  were  generated  once  in  step  2)  is  neces- 
sary to  eliminate  interpolation  errors.   This  feature  however  will  not 
introduce  significantly  more  computations  when  the  control  lies  always 
on  a  constraint. 
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